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Abstract 

These lectures were given at TASI 2012 and are directed at a level suitable for graduate 
students in High Energy Physics. They are intended to give an introduction to the the- 
ory and phenomenology of quantum chromodynamics (QCD), focusing on collider physics 
applications. The aim is to bring the reader to a level where informed decisions can be 
made concerning different approaches and their uncertainties. The material is divided 
into five main areas: 1) fundamentals, 2) fixed-order perturbative QCD, 3) Monte Carlo 
event generators and parton showers, 4) Matching at Leading and Next-to-Leading Order, 
and 5) Soft QCD physics. 
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1 Introduction 

when probed at very short wavelengths, QCD is essentially a theory of free 'partons' — quarks 
and gluons — which only scatter off one another through relatively small quantum corrections, 
that can be systematically calculated. At longer wavelengths, of order the size of the proton 
~ Ifm = 10~^^m, however, we see strongly bound towers of hadron resonances emerge, with 
string-like potentials building up if we try to separate their partonic constituents. Due to our 
inability to perform analytic calculations in strongly coupled field theories, QCD is therefore 
still only partially solved. Nonetheless, all its features, across all distance scales, are believed 
to be encoded in a single one-line formula of alluring simplicity; the Lagrangian of QCD. 

The consequence for collider physics is that some parts of QCD can be calculated in terms 
of the fundamental parameters of the Lagrangian, whereas others must be expressed through 
models or functions whose effective parameters are not a priori calculable but which can be 
constrained by fits to data. However, even in the absence of a perturbative expansion, there 
are still several strong theorems which hold, and which can be used to give relations between 
seemingly different processes. (This is, e.g., the reason it makes sense to constrain parton 
distribution functions in ep collisions and then re-use the same ones for pp collisions.) Thus, 
in the chapters dealing with phenomenological models we shall emphasize that the loss of a 
factorized perturbative expansion is not equivalent to a total loss of predictivity. 

An alternative approach would be to give up on calculating QCD altogether and use leptons 
instead. Formally, this amounts to summing inclusively over strong-interaction phenomena, 
when such are present. While such a strategy might succeed in replacing what we do know 
about QCD by "unit/', however, even the most adamant chromophobe must acknowledge a 
few basic facts of collider physics for the next decade(s): 1) At the Tevatron and LHC, the 
initial states are hadrons, and hence, at the very least, well-understood and precise parton 
distribution functions (PDFs) will be required; 2) high precision will mandate calculations to 
higher orders in perturbation theory, which in turn will involve more QCD; 3) the requirement 
of lepton isolation makes the very definition of a lepton depend implicitly on QCD, and 4) the 
rate of jets that are misreconstructed as leptons in the experiment depends explicitly on it. 
Finally, 5) though many new-physics signals do give observable signals in the lepton sector, 
this is far from guaranteed. It would therefore be unwise not to attempt to solve QCD to the 
best of our ability, the better to prepare ourselves for both the largest possible discovery reach 
and the highest attainable subsequent precision. 

In addition, or perhaps as a consequence, the field of QCD is currently experiencing some- 
thing of a revolution. On the perturbative side, new methods to compute scattering am- 
plitudes with very high particle multiplicities are being developed, together with advanced 
techniques for combining such amplitudes with all-orders resummation frameworks. On the 
non-perturbative side, the wealth of data on soft-physics processes from the LHC is forcing us 
to reconsider the reliability of the standard fragmentation models, and heavy-ion collisions are 
providing new insights into the collective behavior of hadronic matter. The study of cosmic 
rays impinging on the Earth's atmosphere challenges our ability to extrapolate fragmentation 
models from collider energy scales to the region of ultra-high energy cosmic rays. And finally, 
dark-matter annihilation processes in space can produce colored particles either directly or via 
decays of unstable resonances, making predictions of their spectra sensitive to the accuracy of 
the fragmentation modeling. 
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Meson-Nucleon Scattering and Nucleon Isobars* 



Keith A. Brueckner 
Department of Physics, Indiana University, Bloomington, Indiana 
(Received December 17, 1951) 



"[...] It is concluded that the apparently anomalous features of the scattering can be 
interpreted to be an indication of a resonant meson-nucleon interaction corresponding to 
a nucleon isobar with spin |, isotopic spin |, and with an excitation energy of 277 MeV." 

Figure 1: The title and part of the abstract of the 1951 paper [1] (published in 1952) in which 
the discovery of the A++ baryon was announced. 



In the following, we shall focus squarely on QCD for mainstream collider physics. This 
includes factorization, hard processes, infrared safety, parton showers and matching, event 
generators, hadronization, and the so-called underlying event. While not covering everything, 
hopefully these topics can also serve at least as stepping stones to more specialized issues that 
have been left out, such as twistor-inspired techniques, heavy flavors, or forward physics, or 
to topics more tangential to other fields, such as lattice QCD or heavy-ion physics. 

1.1 A First Hint of Color 

Looking for new physics, as we do now at the LHC, it is instructive to consider the story of the 
discovery of color. The first hint was arguably the A++ baryon, found in 1951 [1]. The title 
and part of the abstract from this historical paper are reproduced in figure 1 . In the context of 
the quark model — which first had to be developed, successively joining together the notions 
of spin, isospin, strangeness, and the eightfold way^ — the flavor and spin content of the A++ 



clearly a highly symmetric configuration. However, since the A++ is a fermion, it must have 
an overall antisymmetric wave function. In 1965, fourteen years after its discovery, this was 
finally understood by the introduction of color as a new quantum number associated with the 
group SU(3) [2, 3]. The A++ wave function can now be made antisymmetric by arranging its 
three quarks antisymmetrically in this new degree of freedom. 



hence solving the mystery. 

More direct experimental tests of the number of colors were provided first by measure- 
ments of the decay width of vr*^ —?- 77 decays, which is proportional to N^, and later by the 
famous "R" ratio in e+e~ collisions (i? = a{e~^e~ — )• qq)/a{e~^e~ — )• fi^fi^)), which is pro- 
portional to Nc, see e.g. [4]. Below, in Section 1.2 we shall see how to calculate such color 
factors. 

^In physics, the "eightfold way" refers to the classification of the lowest-lying pseudoscalar mesons and spin-1/2 
baryons within octets in SU(3)-flavor space («,d, s). The is part of a spin- 3/2 baryon decuplet, a "tenfold 
way" in this terminology. 



baryon is: 




(1) 




(2) 
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1.2 The Lagrangian of QCD 

Quantum Chromodynamics is based on the gauge group SU(3), the Special Unitary group in 3 
(complex) dimensions. In the context of QCD, we represent this group as a set of unitary 3x3 
matrices with determinant one. This is called the adjoint representation and can be used to rep- 
resent gluons in color space. Since there are 9 linearly independent unitary complex matrices, 
one of which has determinant —1, there are a total of 8 independent directions in the adjoint 
color space, i.e., the gluons are octets. In QCD, these matrices can operate both on each other 
(gluon self- interactions) and on a set of complex 3-vectors (the fundamental representation), 
the latter of which represent quarks in color space. The fundamental representation has one 
linearly independent basis vector per degree of SU(3), and hence the quarks are triplets. 
The Lagrangian of QCD is 

C = i^iii^niD.hji^i, - m.i^iip,, - \f^,F-^'' , (3) 

where denotes a quark field with color index i, ipq = {iJ-'qR, i^',,c,i^qBf , 7^ is a Dirac matrix 
that expresses the vector nature of the strong interaction, with ^ being a Lorentz vector index, 
niq allows for the possibility of non-zero quark masses (induced by the standard Higgs mech- 
anism or similar), F^^ is the gluon field strength tensor for a gluon with color index a (in the 
adjoint representation, i.e., a e [1, • . . , 8]), and is the covariant derivative in QCD, 

{D^)ij = 6^jd^ - igst%A1 , (4) 

with gs the strong coupling (related to a.^ by = Airas', we return to the strong coupling 
in more detail below), the gluon field with (adjoint-representation) color index a, and 
proportional to the hermitean and traceless Gell-Mann matrices of SU(3), 
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(5) 



These generators are just the SU(3) analogs of the Pauli matrices in SU(2). By convention, the 
constant of proportionality is normally taken to be 



(6) 



This choice in turn determines the normalization of the coupling gs, via equation (4), and fixes 
the values of the SU(3) Casimirs and structure constants, to which we return below. 

An example of the color flow for a quark-gluon interaction in color space is given in fig- 
ure 2. 



1.3 Color Factors 

Typically, we do not measure color in the final state — instead we average over all possible in- 
coming colors and sum over all possible outgoing ones, wherefore QCD scattering amplitudes 
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Figure 2: Illustration of a qqg vertex in QCD, before summing/averaging over colors: a gluon 
in a state represented by A^ interacts with quarks in the states ijjqn and V'gC- 



(squared) in practice always contain sums over quark fields contracted with Gell-Mann matri- 
ces. These contractions in turn produce traces which yield the color factors that are associated 
to each QCD process, and which basically count the number of "paths through color space" 
that the process at hand can take^. 

A very simple example of a color factor is given by the decay process Z — > qq. This vertex 
contains a simple 6ij in color space; the outgoing quark and antiquark must have identical 
(anti-) colours. Squaring the corresponding matrix element and summing over final-state col- 
ors yields a color factor of 

e+e- ^ Z ^qq : ^ \M\'^ cx 5ij5*i = Tr{6} = Nc = S , (7) 

colors 

since i and j are quark (i.e., 3-dimensional fundamental-representation) indices. 

A next-to-simplest example is given hy qq —?- ^* jZ —?■ i'^£~ (usually referred to as the 
Drell-Yan process [5]), which is just a crossing of the previous one. By crossing symmetry, the 
squared matrix element, including the color factor, is exactly the same as before, but since the 
quarks are here incoming, we must average rather than sum over their colors, leading to 

qq^Z^ e+e- : ^ J] \M\^ « ^d,,6*., = ^Tr{J} = ^ , (8) 

colors 

where the color factor now expresses a suppression which can be interpreted as due to the fact 
that only quarks of matching colors are able to collide and produce a Z boson. The chance 
that a quark and an antiquark picked at random from the colliding hadrons have matching 
colors is 1/Nc- 

Similarly, £q — > iq via i-channel photon exchange (usually called Deep Inelastic Scattering 
— DIS — with "deep" referring to a large virtuality of the exchanged photon), constitutes yet 
another crossing of the same basic process, see figure 3. The color factor in this case comes 
out as unity. 

^The convention choice represented by equation (6) introduces a "spurious" factor of 2 for each power of the 
coupHng as. Although one could in principle absorb that factor into a redefinition of the coupling, effectively 
redefining the normalization of "unit color charge", the standard definition of Os is now so entrenched that alter- 
native choices would be counter-productive, at least in the context of a supposedly pedagogical review. 
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Figure 3: Illustration of the three crossings of the interaction of a lepton current (black) with a 
quark current (red) via an intermediate photon or Z boson, with corresponding color factors. 



To illustrate what happens when we insert (and sum over) quark-gluon vertices, such as the 
one depicted in figure 2, we take the process Z — > 3jets. The color factor for this process can 
be computed as follows, with the accompanying illustration showing a corresponding diagram 
(squared) with explicit color-space indices on each vertex: 

Z qgq : 

J2 \M\' oc 5,,t%{t%5:,r 

= :^Tr{<5} = 4, 



colors 





(9) 



where the last Tr{5} = 8, since the trace runs over indices in the 8-dimensional adjoint repre- 
sentation. 

The tedious task of taking traces over SU(3) matrices can be greatly alleviated by use of the 
relations given in Table 1. In the standard normalization convention for the SU(3) generators, 
equation (6), the Casimirs of SU(3) appearing in Table 1 are^ 



Tr 



1 



Cf 



Ca = Nc = 3 



(10) 



In addition, the gluon self-coupling on the third line in Table 1 involves factors of f^^^. These 
are called the structure constants of QCD and they enter due to the non-Abelian term in the 
gluon field strength tensor appearing in equation (3), 



(11) 



'See, e.g., [7, Appendix A.3] for how to obtain the Casimirs in other normalization conventions. 
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Trace Relation 



Indices 



Occurs in Diagram Squared 



ab 



a,6e [1,...,! 



ae [1,...,8] 
z, j,fc G [1, ... ,3] 



a, 6, c, d e [1, ■ ■ ■ , 8] 



(^Sjk^ie- T^SijSke^ ij,k,ie [1,...,3] 




(Fierz) 



Table 1: Trace relations for t matrices (convention-independent). More relations can be found 
in [6, Section 1.2] and in [7, Appendix A.3]. 



Structure Constants of SU(3) 

/l23 = 1 



/l47 = /; 



246 



/257 — /; 



1 



The structure constants of SU(3) are listed 
in the table to the right. Expanding the 
F^pF^" term of the Lagrangian using equa- 
tion (11), we see that there is a 3-gluon and 
a 4-gluon vertex that involve /"''^, the lat- 
ter of which has two powers of / and two 
powers of the coupling. 
Finally, the last line of Table 1 is not re- 
ally a trace relation but instead a useful 
so-called Fierz transformation. It is often 
used, for instance, in shower Monte Carlo 
applications, to assist in mapping between 
color flows in Nc = 3, in which cross sec- 
tions and splitting probabilities are calcu- 
lated, and those in Nc — ^ oo, used to rep- 
resent color flow in the MC "event record". 

A gluon self-interaction vertex is illustrated in figure 4, to be compared with the quark- 
gluon one in figure 2. We remind the reader that gauge boson self-interactions are a hallmark 
of non-Abelian theories and that their presence leads to some of the main differences between 
QED and QCD. One should also keep in mind that the color factor for the vertex in figure 4, 
Ca, is roughly twice as large as that for a quark, Cp- 



/l56 — /367 
/458 = /678 



345 - 
1 

2 



Antisymmetric in all indices 
All other f,jk = 



(12) 
(13) 
(14) 

(15) 



1.4 The Strong Coupling 

To first approximation, QCD is scale invariant. That is, if one "zooms in" on a QCD jet, one will 
find a repeated self-similar pattern of jets within jets within jets, reminiscent of fractals such 
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Figure 4: Illustration of a ggg vertex in QCD, before summing/averaging over colors: interac- 
tion between gluons in the states A^, A^, and A^ is represented by the structure constant /^^^. 



as the famous Mandelbrot set in mathematics, or the formation of frost crystals in physics. 
In the context of QCD, this property was originally called light-cone scaling, or Bjorken scal- 
ing after the physicist James D. Bjorken. This type of scaling is closely related to the class of 
angle-preserving symmetries, called conformal symmetries. In physics today, the terms "confor- 
mal" and "scale invariant" are used interchangeably'^. Conformal invariance is a mathematical 
property of several QCD-"like" theories which are now being studied (such as = 4 super- 
S3mimetric relatives of QCD). It is also closely related to the physics of so-called "unparticles", 
though that is a relation that goes beyond the scope of these lectures. 

Regardless of the labeling, if the strong coupling did not run (we shall return to the running 
of the coupling below), Bjorken scaling would be absolutely true. QCD would be a theory with 
a fixed coupling, the same at all scales. This simplified picture already captures some of the 
most important properties of QCD, as we shall discuss presently. 

In the limit of exact Bjorken scaling — QCD at fixed coupling — properties of high-energy 
interactions are determined only by dimensionless kinematic quantities, such as scattering 
angles (pseudorapidities) and ratios of energy scales^. For applications of QCD to high- 
energy collider physics, an important consequence of Bjorken scaling is thus that the rate 
of bremsstrahlung jets, with a given transverse momentum, scales in direct proportion to 
the hardness of the fundamental partonic scattering process they are produced in association 
with. For instance, in the limit of exact scaling, a measurement of the rate of 10-GeV jets 
produced in association with an ordinary Z boson could be used as a direct prediction of the 
rate of 100-GeV jets that would be produced in association with a 900-GeV Z' boson, and so 
forth. Our intuition about how many bremsstrahlung jets a given type of process is likely to 
have should therefore be governed first and foremost by the ratios of scales that appear in that 
particular process, as has been highlighted in a number of studies focusing on the mass and pj_ 
scales appearing, e.g., in Beyond-the- Standard-Model (BSM) physics processes [8, 9, 10, 11]. 
Bjorken scaling is also fundamental to the understanding of jet substructure in QCD, see, e.g., 

"^Strictly speaking, conformal symmetry is more restrictive than just scale invariance, but examples of scale- 
invariant field theories that are not conformal are rare. 

^Originally, the observed approximate agreement with this was used as a powerful argument for pointlike 
substructure in hadrons; since measurements at different energies are sensitive to different resolution scales, 
independence of the absolute energy scale is indicative of the absence of other fundamental scales in the problem 
and hence of pointlike constituents. 
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[12, 13]. 

In real QCD, the coupling runs logarithmically with the energy, 

where the function driving the energy dependence, the heta function, is defined as 

/3(as) = -al{bQ + biUg + 6205 + ■■■), (17) 
with LO (1-loop) and NLO (2-loop) coefficients 

, _ llCA-^TRUf 



(18) 



UCl-lOTRCAnf-GTRCpnf 153 -19n^ ^^^^ 
= 24^^ = 24vr^ ' ^'^^ 

In the bo coefficient, the first term is due to gluon loops while the second is due to quark ones. 
Similarly, the first term of the bi coefficient arises from double gluon loops, while the second 
and third represent mixed quark-gluon ones. At higher loop orders, the bi coefficients depend 
explicitly on the renormalization scheme that is used. A brief discussion can be found in the 
PDG review on QCD [14], with more elaborate ones contained in [4, 6]. Note that, if there are 
additional colored particles beyond the Standard-Model ones, loops involving those particles 
enter at energy scales above the masses of the new particles, thus modifying the running of 
the coupling at high scales. This is discussed, e.g., for supersymmetric models in [15]. 

Numerically, the value of the strong coupling is usually specified by giving its value at the 
specific reference scale = M|, from which we can obtain its value at any other scale by 
solving equation (16), 

a.(Q2) = a,(M|) -^--^ , (20) 

l + 6o«.(M|)ln^ + 0(a2) 

with relations including the 0{al) terms available, e.g., in [6]. Relations between scales 
not involving M| can obviously be obtained by just replacing M| by some other scale Q'^ 
ever5where in equation (20). A comparison of running at one- and two-loop order, in both 
cases starting from as{Mz) = 0.12, is given in figure 5. As is evident from the figure, the 
2-loop running is somewhat faster than the 1-loop one. 

As an application, let us prove that the logarithmic running of the coupling implies that 
an intrinsically multi-scale problem can be converted to a single-scale one, up to corrections 
suppressed by two powers of Og, by taking the geometric mean of the scales involved. This 
follows from expanding an arbitrary product of individual factors around an arbitrary scale 
fi, using equation (20), 



i=i ^ ' 

dl{ix) U + 60 a. In f , f " ^ + 0{al) ) , (21) 
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Figure 5: Illustration of the running of at 1- (open circles) and 2-loop order (filled circles), 
starting from the same value of as{Mz) = 0.12. 

whereby the specific single-scale choice /u" = /ii/i2 • • • /Un (the geometric mean) can be seen to 
push the difference between the two sides of the equation one order higher than would be the 
case for any other combination of scales^. 

The appearance of the number of flavors, nj, in implies that the slope of the running 
depends on the number of contributing flavors. Since full QCD is best approximated by n/ = 3 
below the charm threshold, by = 4 from there to the h threshold, and by n/ = 5 above 
that, it is therefore important to be aware that the running changes slope across quark flavor 
thresholds. Likewise, it would change across the threshold for top or for any colored new- 
physics particles that might exist, with a magnitude depending on the particles' color and spin 
quantum numbers. 

The negative overall sign of equation (17), combined with the fact that 6o > (for nj < 
16), leads to the famous result^ that the QCD coupling effectively decreases with energy, called 
asymptotic freedom, for the discovery of which the Nobel prize in physics was awarded to 
D. Gross, H. Politzer, and F. Wilczek in 2004. An extract of the prize announcement runs as 
follows: 

^In a fixed-order calculation, the individual scales jii, would correspond, e.g., to the n hardest scales appearing 
in an infrared safe sequential clustering algorithm applied to the given momentum configuration. 

^ Perhaps the highest pinnacle of fame for equation (17) was reached when the sign of it featured in an episode 
of the TV series "Big Bang Theory". 
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What this year's Laureates discovered was something that, at first sight, seemed 
completely contradictory. The interpretation of their mathematical result was 
that the closer the quarks are to each other, the weaker is the "color charge". 
When the quarks are really close to each other, the force is so weak that they 
behave almost as free particle^ . This phenomenon is called "asymptotic free- 
dom". The converse is true when the quarks move apart: the force becomes 
stronger when the distance increased. 



"More correctly, it is the coupling rather than the force which becomes weak as the distance 
decreases. The 1 /r^ Coulomb singularity of the force is only dampened, not removed, by the 
diminishing coupling. 

'More correctly, it is the potential which grows, linearly, while the force becomes constant. 

Among the consequences of asymptotic freedom is that perturbation theory becomes better 
behaved at higher absolute energies, due to the effectively decreasing coupling. Perturbative 
calculations for our 900-GeV Z' boson from before should therefore be slightly faster con- 
verging than equivalent calculations for the 90-GeV one. Furthermore, since the running of 
as explicitly breaks Bjorken scaling, we also expect to see small changes in jet shapes and in 
jet production ratios as we vary the energy. For instance, since high-p^ jets start out with 
a smaller effective coupling, their intrinsic shape (irrespective of boost effects) is somewhat 
narrower than for low-pi^ jets, an issue which can be important for jet calibration. Our current 
understanding of the running of the QCD coupling is summarized by the plot in figure 6, taken 
from a recent comprehensive review by S. Bethke [16]. 

As a final remark on asymptotic freedom, note that the decreasing value of the strong 
coupling with energy must eventually cause it to become comparable to the electromagnetic 
and weak ones, at some energy scale. Beyond that point, which may lie at energies of order 
10^^ — 10^^ GeV (though it may be lower if as yet undiscovered particles generate large correc- 
tions to the running), we do not know what the further evolution of the combined theory will 
actually look like, or whether it will continue to exhibit asymptotic freedom. 

Now consider what happens when we run the coupling in the other direction, towards 
smaller energies. Taken at face value, the numerical value of the coupling diverges rapidly at 
scales below 1 GeV, as illustrated by the curves disappearing off the left-hand edge of the plot 
in figure 6. To make this divergence explicit, one can rewrite equation (20) in the following 
form, 

= — ^ , (22) 
bo In^ 

where 

A~ 200MeV (23) 

specifies the energy scale at which the perturbative coupling would nominally become infinite, 
called the Landau pole. (Note, however, that this only parametrizes the purely perturbative 
result, which is not reliable at strong coupling, so equation (22) should not be taken to imply 
that the physical behavior of full QCD should exhibit a divergence for Q — > A.) 

Finally, one should be aware that there is a multitude of different ways of defining both A 
and as{Mz)- At the very least, the numerical value one obtains depends both on the renor- 
malization scheme used (with the dimensional-regularization-based "modified minimal sub- 
traction" scheme, MS, being the most common one) and on the perturbative order of the 
calculations used to extract them. As a rule of thumb, fits to experimental data typically yield 
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Figure 6: Illustration of the running of as in a theoretical calculation (yellow shaded band) 
and in physical processes at different characteristic scales, from [16]. 



smaller values for as(Mz) the higher the order of the calculation used to extract it (see, e.g., 
[16, 17]), with as{Mz)\LO Z as{Mz)\NLO Z as{Mz)\N]<SLO- Further, since the number of 
flavors changes the slope of the running, the location of the Landau pole for fixed as{Mz) 
depends explicitly on the number of flavors used in the running. Thus each value of n/ is 
associated with its own value of A, with the following matching relations across thresholds 
guaranteeing continuity of the coupling at one loop, 

_2_ _2_ 

n^ = 4o5 : A5 = A4 " A4 = A5 " , (24) 

(Ao \ 25 / rr? \ 27 

^) ^'^^^[li) ■ ^^^^ 

It is sometimes stated that QCD only has a single free parameter, the strong coupling. 
However, even in the perturbative region, the beta function depends explicitly on the number 
of quark flavors, as we have seen, and thereby also on the quark masses. Furthermore, in 
the non-perturbative region around or below Aqcd, the value of the perturbative coupling, as 
obtained, e.g., from equation (22), gives little or no insight into the behavior of the full theory. 
Instead, universal functions (such as parton densities, form factors, fragmentation functions, 
etc), effective theories (such as the Operator Product Expansion, Chiral Perturbation Theory, 
or Heavy Quark Effective Theory), or phenomenological models (such as Regge Theory or the 
String and Cluster Hadronization Models) must be used, which in turn depend on additional 
non-perturbative parameters whose relation to, e.g., as{Mz), is not a priori known. 

For some of these questions, such as hadron masses, lattice QCD can furnish important 
additional insight, but for multi-scale and/or time-evolution problems, the applicability of 
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lattice methods is still severely restricted; the lattice formulation of QCD requires a Wick 
rotation to Euclidean space. The time-coordinate can then be treated on an equal footing with 
the other dimensions, but intrinsically Minkowskian problems, such as the time evolution of a 
system, are inaccessible. The limited size of current lattices also severely constrain the scale 
hierarchies that it is possible to "fit" between the lattice spacing and the lattice size. 



— 14 — 



p. Skands 



Introduction to QCD 




Figure 7: Left: Rutherford scattering of quarks in QCD, exemplifying the type of process that 
dominates the short-distance interaction cross section at hadron colliders. Right: an example 
of what such a reaction looks like in a detector, in this case the ATLAS experiment. 



2 Hard Processes 

Our main tool for solving QCD at high energy scales, Q » Aqcd, is perturbative quantum 
field theory, the starting point for which is Matrix Elements (MEs) which can be calculated 
systematically at fixed orders (FO) in the strong coupling ag. At least at lowest order (LO), the 
procedure is standard textbook material [7] and it has also by now been highly automated, by 
the advent of tools like Madgraph [18], Calchep [19] / Comphep [20], and several others 
[21, 22, 23, 24, 25, 26, 27]. Here, we require only that the reader has a basic familiarity with 
the methods involved from graduate-level particle physics courses based, e.g., on [7, 4]. Our 
main concern are the uses to which these calculations are put, their limitations, and ways to 
improve on the results obtained with them. 

For illustration, take one of the most commonly occurring processes in hadron collisions: 
Rutherford scattering of two quarks via a t-channel gluon exchange — figure 7 — which at 
leading order has the differential cross section 



, , dcr _ vr 4 2 ■s^ + 



qq ^ qq : —r = ^7:as — ^ — , (26) 



with the 2 — ;> 2 Mandelstam variables ("hatted" to emphasize that they refer to a partonic 
2 — ;> 2 scattering rather than the full pp — )- jets process) 

s = {Pi+P2f, (27) 

n9 ^ (1 — COS 9) 
t = {P3-Pi? = -s^- ^ (28) 

f \2 .(l + cos6') 
u = {pA-piY = -s^ ^ . (29) 

Reality, however, is more complicated; the picture on the right-hand pane of figure 7 shows 
a real dijet event, as recorded by the ATLAS experiment. The complications to be addressed 
when going from left to right in figure 7 are: firstly, additional jets, a.k.a. real-emission cor- 
rections, which can significantly change the topology of the final state, potentially shifting jets 
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Figure 8: Illustration of partonic fluctuations inside a proton beam (from [28]). 



in or out of an experimentally defined acceptance region. Secondly, loop factors, a.k.a. virtual 
corrections, change the number of available quantum paths through phase space, and hence 
modify the normalization of the cross section (total and differential). And finally, additional 
corrections are generated by confinement and by the so-called underlying event. These cor- 
rections must be taken into account to complete our understanding of QCD and connect the 
short-distance physics with macroscopic experiments. Apart from the perturbative expansion 
itself, the most powerful tool we have to organize this vast calculation, is factorization. 



2. 1 Factorization 



In high-energy scattering problems involving hadrons in the initial state, we immediately face 
the complication that hadrons are composite, with a time-dependent structure illustrated in 
figure 8; there are partons within clouds of further partons, constantly being emitted and ab- 
sorbed. Thus, before we can use perturbatively calculated partonic scattering matrix elements, 
we must first address the partonic structure of the colliding hadron(s). 

For the hadron to remain intact, the fluctuations inside it must involve momentum transfers 
smaller than the confinement scale. Indeed, high-virtuality fluctuations are suppressed by 
powers of 




(30) 



with A the confinement scale (~ 200 MeV, see section 1.4) and |A;| the virtuality of the fluctu- 
ation. Thus, most fluctuations occur over timescales ~ 1/A. 

A hard perturbative probe, on the other hand, such as the exchanged photon in DIS (fig- 
ure 3), interacts over a much shorter timescale 1/Q <C 1/A, during which the partonic fluctu- 
ations in the struck hadron appear almost frozen. The hard probe effectively takes an instan- 
taneous snapshot of the hadron structure, at a characteristic resolution given by ~ 1/Q. 

This is formalized by the factorization theorem [29] (see also the TASI lectures by George 
Sterman [30]), which expresses the independence of long-wavelength (soft) structure on the 
nature of the hard (short-distance) process. Originally formulated for DIS, factorization allows 
us to write the cross section for lepton-hadron scattering as a convolution of a non-perturbative 
but universal (i.e., process-independent) parton density function (PDF) and a perturbatively 
calculable partonic scattering cross section. Denoting the fraction of the hadron momentum 
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carried by parton i by Xi 



Pi = XiPh 



(31) 



we may write the lepton-hadron cross section on factorized form (see, e.g., [31, 4]), 



with i an index running over all possible parton types^ in the incoming hadron and / enumer- 
ating all possible (partonic) final states, with Lorentz-invariant phase space, <I>y. 

The parton density functions (PDFs), fi/^, parametrize the distribution of partons inside the 
target hadron. They are not a priori calculable and must be constrained by fits to data. This is 
discussed in section 2.2. 

The partonic cross section, da, knows nothing of the target hadron apart from the fact that 
it contained the struck parton. It is calculable within perturbation theory, as will be discussed 
in section 2.3. 

The dividing line between the two is drawn at an arbitrary ("user-defined") scale fip, called 
the factorization scale. There is some arbitrariness involved in choosing a value for fip. Some 
heuristic arguments to guide in the choice of factorization scale are the following. On the 
long-distance side, the PDFs include a (re) summation of fluctuations inside fluctuations up 
to virtualities of order /xi?. It would therefore not make much sense to take /ip significantly 
larger than the scales characterizing resolved particles on the short-distance side of the calcu- 
lation (i.e., the particles appearing explicitly in <I>^); otherwise the PDFs would be including 
sums over fluctuations that happen on timescales shorter than those probed by the physical 
process. Similarly, fip should also not be taken much lower than the scale (s) appearing in the 
hard process. For matrix elements characterized by a single well-defined scale, such as the 

scale in DIS or the invariant-mass scale s in Drell-Yan production (qq — > Z/j* — )- £^£~), 
such arguments essentially fix the preferred scale choice to fip = Q or fip = ^/§, respectively, 
which may then be varied by a factor of 2 (or larger) around the nominal value in order to 
estimate uncertainties. For multi-scale problems, however, such as pp — ;> Z/W + njets, there 
are several a priori equally good choices available, from the lowest to the highest QCD scales 
that can be constructed from the final-state momenta, usually with several dissenting groups 
of theorists arguing over which particular choice is best. Suggesting that one might simply 
measure the scale would not really be an improvement, as the factorization scale is fundamen- 
tally unphysical and therefore unobservable (similarly to gauge or convention choices). One 
plausible strategy is to look at higher-order (NLO or NNLO) calculations, in which correction 
terms appear that cancel the dependence on the scale choice, stabilizing the final result. From 
such comparisons, a "most stable" initial scale choice can in principle be determined, which 
then furnishes a reasonable starting point, but we emphasize that the question is intrinsically 
ambiguous, and no golden recipe is likely to magically give all the right answers. The best we 
can do is to vary the value of //i? not only by an overall factor, but also by exploring different 
possible forms for its functional dependence on the momenta appearing in <I> ^ . A complemen- 
tary useful discussion of the pros and cons of different factorization scale choices can be found 
in the TASI lectures by Tilman Plehn [32]. 

^Typically, only quarks and gluons are included in this sum, but also photons and even leptons can in principle 
be included. Similarly, parton density functions are normally used to describe hadrons, but can also be defined, 
e.g., to describe the cloud of virtual photons (and fermion pairs) surrounding an electron. 




(32) 
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Secondly, and more technically, at NLO and beyond one also has to settle on a factorization 
scheme in which to do the calculations. For all practical purposes, students focusing on LHC 
physics are only likely to encounter one such scheme, the modified minimal subtraction (MS) 
one already mentioned in the discussion of the definition of the strong coupling in Section 1.4. 
At the level of these lectures, we shall therefore not elaborate further on this choice here. 

We note that factorization can also be applied multiple times, to break up a complicated 
calculation into simpler pieces that can be treated as approximately independent. This will be 
very useful when dealing with successive emissions in a parton shower, section 3.2, or when 
factoring off decays of long-lived particles from a hard production process, section 3.4. 

We round off the discussion of factorization by mentioning a few caveats the reader should 
be aware of. (See [30] for a more technical treatment.) 

Firstly, the proof only applies to the first term in an operator product expansion in "twist" 
= mass dimension - spin. Since operators with higher mass dimensions are suppressed by the 
hard scale to some power, this leading twist approximation becomes exact in the limit Q — ^ oo, 
while at finite Q it neglects corrections of order 

[ln(Q7A2)]'™<2" 

In section 5, we shall discuss some corrections that go beyond this approximation, in the 
context of multiple parton-parton interactions. 

Secondly, the proof only really applies to inclusive cross sections in DIS [29] and in Drell- 
Yan [33]. For all other hadron-initiated processes, factorization is an ansatz. For a general 
hadron-hadron process, we write the assumed factorizable cross section as: 

df^/^l/l2 = / Yl / fi/i^^ fi/h2 i^j^l^p) . ^r^L ■ (34) 

• ■ r J ^ .1 t 



Higher Twist : ^ 'J^ (n = 2 for DIS) . (33) 



Note that, if da is divergent (as, e.g., Rutherford scattering is) then the integral over d^^ 
must be regulated, e.g. by imposing some explicit minimal transverse-momentum cut and/ or 
other phase-space restrictions. 



2.2 Parton Densities 

The parton density function, fi/h{x,iJ-'j?), represents the effective density of partons of type/flavor 
i, as a function of the momentum fraction^ Xi, when a hadron of t5^e h is probed at the fac- 
torization scale fxp- The PDFs are non-perturbative functions which are not a priori calculable, 
but a perturbative differential equation governing their evolution with fxp can be obtained 
by requiring that physical scattering cross sections, such as the one for DIS in equation (32), 
be independent of fip to the calculated orders [34]. The resulting renormalization group 
equation (RGE) is called the DGLAP^° equation and can be used to "run" the PDFs from one 
perturbative resolution scale to another (its evolution kernels are the same as those used in 
parton showers, to which we return in section 3.2). 

This means that we only need to determine the form of the PDF as a function of x a single 
(arbitrary) scale, fiQ- We can then get its form at any other scale fip^Y simple RGE evolution. 

'Recall: the x fraction is defined in equation (31). 
^"DGLAP: Dokshitzer-Gribov-Lipatov-Altarelli-Parisi. 



— 18 — 



p. Skands 



Introduction to QCD 



0**2= 4 GeVtt2 

up MSTW2008L0 

upbar MSTW2008L0 

chorm MSTW2008L0 

.... gljon MSTW2008L0 x 0. 




Qt,2= 10000 GeV«*2 

up MSTW2008LO 

upbar MSTW2008LO 

charm USTW2008LO 

.... gluon MSTW2008LO x 0.1 




Figure 9: Illustration of the change of the u (black), u (red, dashed), c (green, dotted), and 
g (blue, dot-dashed) distributions, from Q = = 2GeV (left) to Q = = 100 GeV 
(right). Note that a factor 0.1 has been applied to the gluon distribution. Plots made using the 
HEPDATA online tool [40]. 



In the context of PDF fits (constraining the form of the PDF functions by fitting cross sections 
to experimental data, e.g., from DIS [35, 36], Drell-Yan [37, 38], and pp — ^ jets [39]), the 
reference scale fiQ is usually taken to be relatively low, of order one or a few GeV. 

The behavior of the PDFs as we evolve ^ip from a low scale, 2 GeV, to a high one, 100 GeV, 
is illustrated in figure 9, for the MSTW^i 2008 LO^^ p^p ggt [41]. At low Q = ^u^ = 2 GeV 
(left), the proton structure is dominated by a few hard quarks (a "valence bump" is clearly 
visible around x ~ 0.2), while at higher scales Q = 100 GeV (right) we predominantly resolve 
fluctuations within fluctuations, yielding increasingly large gluon- and sea-quark distributions 
with rather small x values, while the valence quarks play a progressively smaller role. 

We note that different collaborations, like CTEQ, MSTW, NNPDF, etc., use different ansatze 
for the form of /(x, //g)- They may also include different data in the fits, and/ or treat or weight 
the data differently. Thus, results from different groups may not always be mutually compati- 
ble. An example is given in figure 10, which shows the difference between the CTEQ6L1 gluon 
PDF [42] (red dashed) and the MSTW 2008 LO one [41], normalized to MSTW (which would 
thus be a flat line at zero), at /xi? = 10 GeV. The y axis shows the relative difference between 
the sets, in per cent. Also shown are the 90% CL contours computed from the uncertainty vari- 
ations included in the MSTW 2008 LO set (black). Using only the MSTW uncertainty band, 
one would arrive at an estimated ~ 5% uncertainty over most of the x range, while including 
the CTEQ6L1 set would increase that to ~ 10%. At NLO, this discrepancy is reduced, but not 
removed. A significant effort is currently being undertaken within the PDF community to agree 
on common, and more comprehensive, ways of defining PDF uncertainty bands [39, 43]. This 

"MSTW: Martin-Stirling-Thorne-Watt. 

^^The "LO" means that the fit was performed using LO matrix elements in the cross section formulae. 
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Figure 10: Illustration of the difference between the MSTW 2008 and CTEQ6 LO gluon PDFs 
at /ii? = lOGeV. All curves are normalized to the central MSTW 2008 prediction. The black 
solid lines show the 90% CL MSTW variations, while the dashed red line shows the CTEQ6L1 
distribution. 

is complicated due to the different ways of defining f{x, fi'^) and due to the experimental data 
sets not always being fully compatible with one another For the time being, it is recommended 
to try at least sets from two different groups, for a comprehensive uncertainty estimate. 

Occasionally, the words structure functions and parton densities are used interchangeably. 
However, there is an important distinction between the two, which we find often in (quantum) 
physics: the former is a physical observable used to parametrize the DIS cross sections (see 
G-g- [4]), while the latter is a "fundamental" quantity extracted from it. In particular, since the 
parton densities are not, themselves, physically observable, they can only be defined within 
a specific factorization scheme, order by order in perturbation theory. The only exception 
is at leading order, at which they have the simple physical interpretation of parton number 
densities. When going to higher orders, we tend to keep the simple intuitive picture from 
LO in mind, but one should be aware that the fundamental relationship between PDFs and 
measured quantities is now more complicated (due to the interplay between the PDFs and the 
real and virtual corrections to the LO cross section), and that the parton densities no longer 
have a clear probabilistic interpretation starting from NLO. 

The reader should also be aware that there is some ambiguity whether NLO PDFs should 
be used for LO calculations. In principle, the higher-order PDFs are better constrained and 
the difference between, e.g., an NLO and an LO set should formally be beyond LO precision, 
so that one might be tempted to simply use the highest-order available PDFs for any calcula- 
tion. However, higher-order terms can sometimes be absorbed, at least partially, into effective 
lower-order coefficients. In the context of PDFs, the fit parameters of lower-order PDFs will 
attempt to compensate for missing higher-order contributions in the matrix elements. To the 
extent those higher-order contributions are universal, this is both desirable and self-consistent. 
This leads to some typical qualitative differences between LO and NLO PDFs, illustrated in 
figure 11: NLO PDFs tend to be smaller at low x and slightly larger at high x, than LO ones. 
Thus, it is quite possible that using an NLO PDF in conjunction with LO matrix elements can 
give a worse agreement with data than LO PDFs do. 
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Figure 11: Illustration of the change between PDF fits using LO and NLO matrix elements: 
the g distribution at LO G^lack) and NLO (red, dashed), and the u distribution at LO (green, 
dotted) and NLO (blue, dot-dashed), for the MSTW 2008 PDF sets [41], at Q = /ii. = 10 GeV. 
Plots made using the HEPDATA online tool [40] . 

Finally, another oft-raised question concerns which PDF sets to use for the parton-shower 
evolution in Monte Carlo generators. Importantly, the equations driving the initial-state show- 
ers in Monte Carlo models are only sensitive to ratios of PDFs [44] . Since the shower evolution 
typically only has leading-logarithmic (LL) precision, it should be theoretically consistent to 
use any (LO or better) PDF set to drive the evolution. However, similarly to above, there will 
be subleading differences between different choices, and one is justified in worrying about 
the level of physical effects that could be generated. Unfortunately, there is currently no way 
to ensure 100% self-consistency. Since PDF fits are not done with MC codes, but instead 
use analytical resummation models (see, e.g., the TASI lectures by Sterman [30]), which are 
not identical to their MC counterparts, the PDF fits are essentially "tuned" to a slightly dif- 
ferent resummation than that incorporated in a given MC model. In practice, not much is 
known about the size and impact of this ambiguity [45]. Known differences include: the size 
of phase space (purely coUinear massless PDF evolution vs. the finite-transverse-momentum 
massive MC phase space), the treatment of momentum conservation and recoil effects, ad- 
ditional higher-order effects explicitly or implicitly included in the MC evolution, choice of 
renormalization scheme and scale, and, for those MC algorithms that do not rely on coUinear 
(DGLAP, see [4]) splitting kernels (e.g., the various kinds of dipole evolution algorithms, see 
[46]), differences in the effective factorization scheme. 

As a baseline, we recommend simply using whatever PDF set the given MC model was 
originally tuned with, since this should de facto (by fitting the available data) reabsorb as 
much of the inconsistency as possible. Furthermore, it should be emphasized that underlying- 
event and minimum-bias models based on multi-parton interactions (see section 5.2) usually 
make the explicit assumption that the PDFs can be interpreted as physical number densities 
even down to very low Q and x, a property which is generally only true for LO PDFs. It must 
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therefore be strongly discouraged to use (N)NLO PDF sets in this context. 
2.3 Fixed-Order QCD 

Consider the production of an arbitrary final state, F (e.g., a Higgs boson, a tt pair, etc). 
Schematically, we may express the (perturbative) all-orders differential cross section for an 
observable O, in the following way: 



where, for compactness, we have suppressed all PDF and luminosity normalization factors. 
■^^F+k amplitude for producing F in association with k additional final-state partons, 

"legs", and with i additional loops. The sums start at = and £ = 0, corresponding to 
the Leading Order for producing F, while higher terms represent real and virtual corrections, 
respectively. 

The purpose of the 6 function is to project out hypersurfaces of constant value of O in the 
full d<I>^_,_^ phase space, with 0{^p^i^) a function that defines O evaluated on each specific 
momentum configuration, (Without the 5 function, the formula would give the total 

integrated cross section, instead of the cross section differentially in C) 

We recover the various fixed-order truncations of perturbative QCD (pQCD) by limiting the 
nested sums in equation (35) to include only specific values of k + £. Thus, 

k = 0, i = Leading Order (usually tree-level) for F production 

k = n, £ = =^ Leading Order for F + n jets 

k + £<n, =^ N"LO for F (includes N"-iLO for + Ijet, N"-2lO for + 



For > 1, we are not considering inclusive F production; instead, we are considering the 
process F + k jets. If we simply integrate over all momenta, as implied by the integration over 
d^F+k equation (35), we would be including configurations in which one or more of the 
k partons are coUinear or soft. Such configurations are infrared divergent in QCD and hence 
must be regulated. Since we talk about collinear and soft divergences (the origins of which will 
be discussed in more detail in sections 2.4 and 3.2), cuts on angles and energies and/or cuts 
on combinations, like transverse momenta, can be used to cut away the problematic regions of 
phase space. 

Recall, however, that pQCD is approximately scale invariant. This implies that a regu- 
larization cut on a dimensionful quantity, like energy or transverse momentum, should be 
formulated as a ratio of scales, rather than as an absolute number For example, a jet with 
p± = 50 GeV would be considered hard and well-separated if produced in association with an 
ordinary Z boson (with hard scale Mz = 91.2 GeV), while the same jet would be considered 
soft if produced in association with a 900-GeV Z' boson (see [8, 9] for more explicit examples). 

The essence of the point is that, if the regularization scale is taken too low, logarithmic 




(35) 



s legs 



s loops 



2jets, and so on up to LO for F + ?ijets) . 
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Figure 12: Coefficients of the perturbative series covered by LO calculations. Left: F pro- 
duction at lowest order Right: F + 2 jets at LO, with the half-shaded box illustrating the 
restriction to the region of phase space with exactly 2 resolved jets. The total power of as for 
each coefficient is n = k + £. (Photo of Max Born from nobelprize . org). 



enhancements of the type 

a-ln™<2"(^^^ (36) 

will generate progressively larger corrections, order by order, which will spoil any fixed-order 
truncation of the perturbative series. Here, Qp is the hard scale associated with the process 
under consideration, while is the scale associated with an additional parton, k. 

A good rule of thumb is that if ak+i ~ o-fc (at whatever order you are calculating), then the 
perturbative series is converging too slowly for a fixed-order truncation of it to be reliable. For 
fixed-order perturbation theory to be applicable, you must place your cuts on the hard process 
such that cTfc+i ^ cjfc. In the discussion of parton showers in Section 3.2, we shall see how the 
region of applicability of perturbation theory can be extended. 

The virtual amplitudes, for £ > I, are divergent for any point in phase space. However, 
as encapsulated by the famous KLN theorem [47, 48], unitarity (which essentially expresses 
probability conservation) puts a powerful constraint on the IR divergences^^, forcing them to 
cancel exactly against those coming from the unresolved real emissions that we had to cut out 
above, order by order, making the complete answer for fixed k+£ = n finite. Nonetheless, since 
this cancellation happens between contributions that formally live in different phase spaces, 
a main aspect of loop-level higher-order calculations is how to arrange for this cancellation 
in practice, either analytically or numerically, with many different methods currently on the 
market. We shall discuss the idea behind subtraction approaches in section 2.4. 

A convenient way of illustrating the terms of the perturbative series that a given matrix- 
element-based calculation includes is given in figure 12. In the left-hand pane, the shaded 
box corresponds to the lowest-order "Born-level" matrix element squared. This coefficient 
is non-singular and hence can be integrated over all of phase space, which we illustrate by 
letting the shaded area fill all of the relevant box. A different kind of leading-order calculation 
is illustrated in the right-hand pane of figure 12, where the shaded box corresponds to the 
lowest-order matrix element squared for F + 2 jets. This coefficient diverges in the part of 
phase space where one or both of the jets are unresolved (i.e., soft or coUinear), and hence 

^^The loop integrals also exhibit UV divergences, but these are dealt with by renormalization. 
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Figure 13: Coefficients of the perturbative series covered by NLO calculations. Left: F produc- 
tion at NLO. Right: F + 1 jet at NLO, with half-shaded boxes illustrating the restriction to the 
region of phase space with exactly 1 resolved jet. The total power of as for each coefficient is 
k + i. 



n 



integrations can only cover the hard part of phase space, which we reflect by only shading the 
upper half of the relevant box. 

Figure 13 illustrates the inclusion of NLO virtual corrections. To prevent confusion, first a 
point on notation: by we intend 

a(^) = / d$o 2Re[Mi'^MS'^*] , (37) 



which is of order relative to the Born level. Compare, e.g., with the expansion of equa- 
tion (35) to order k + i = 1. In particular, o-q^-' should not be confused with the integral over 
the 1-loop matrix element squared (which would be of relative order and hence forms part 

(2) 

of the NNLO coefficient CTq ). Returning to figure 13, the unitary cancellations between real 
and virtual singularities imply that we can now extend the integration of the real correction in 
the left-hand pane over all of phase space, while retaining a finite total cross section, 

J J J (38) 

with (Tq^-* the finite Born-level cross section, and the positive divergence caused by integrating 
the second term over all of phase space is canceled by a negative one coming from the inte- 
gration over loop momenta in the third term. One method for arranging the cancellation of 
singularities — subtraction — is discussed in section 2.4. 

However, if our starting point for the NLO calculation is a process which already has a 
non-zero number of hard jets, we must continue to impose that at least that number of jets 
must still be resolved in the final-state integrations. 



mill "^P_L1 ^P±min '^P.Lmin 



(yf\p± > P±min) + ^(P±l > P±min) + CF^i^ (p± > P±min) , 



(39) 
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Figure 14: Coefficients of the perturbative series covered by an NNLO calculation. The total 
power of as for each coefficient is n = k + i. Green shading represents the full perturbative 
coefficient at the respective k and £. 



where the restriction to at least one jet having p± > p±mm has been illustrated in the right- 
hand pane of figure 13 by shading only the upper part of the relevant boxes. In the second 
term in equation (39), the notation p±i is used to denote that the integral runs over the phase 
space in which at least one "jet" (which may consist of one or two partons) must be resolved 
with respect to p^min- Here, therefore, an explicit dependence on the algorithm used to define 
"a jet" enters for the first time. This is discussed in more detail in the 2009 ESHEP lectures by 
Salam [49]. 

To extend the integration to cover also the case of 2 unresolved jets, we must combine the 
left- and right-hand parts of figure 13 and add the new coefficient 

= |A^«P + 2Re[>l(^)<h, (40) 
as illustrated by the diagram in figure 14. 



2.4 The Subtraction Idea 

According to the KLN theorem, the IR singularities coming from integrating over coUinear and 
soft real-emission configurations should cancel, order by order, by those coming from the IR 
divergent loop integrals. This implies that we should be able to rewrite e.g. the NLO cross 
section, equation (38), as 

(jNLO ^ ^Born ^ Finite | ^ d<^> IM^^^l"^^ + Finite | ^ d^p 2Re[M^p'^ mP*]^^!) 

with the second and third terms having had their common (but opposite-sign) singularities 
canceled out and some explicitly finite quantities remaining. 

The first step towards this goal is to classify all IR singularities that could appear in the 
amplitudes. We know that the IR limits are universal, so they can be classified using a set of 
process-independent functions that only has to be worked out once and for all. A widely used 
such set of functions are the Catani-Seymour (CS) dipole ones [50, 51], a method which by 
now has even been partially automated [52, 53]. Here, we shall instead use a formalism based 
on antennae [54, 55, 56]. The distinction between the two is basically that one antenna is 
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made up of two dipole "ends", hence the antenna formalism tends to generate somewhat fewer 
terms. At NLO, however, there is no fundamental incompatibility — the antennae we use here 
can always be partitioned into two dipole ends, if so desired. (Note: only the antenna method 
has been successfully generalized to NNLO [57, 58]. Other NNLO techniques, not covered 
here, are sector decomposition, see [59, 60], and the generic formalism for hadroproduction of 
colorless states presented in [61].) 

At NLO, the idea with subtraction is thus to rewrite the NLO cross section by adding and 
subtracting a simple function, das' , that encapsulates all the IR limits, 

NLO ^Born , /"j^ f\ XyfC^) |2 ^^NLO 



Finite by Universality 
+ /d..2Re[>,ii)>,i;'-l./d.,,.d.,r . (42) 

^ V ' 

Finite by KLN 

The task now is to construct a suitable form for das ■ A. main requirement is that it should be 
sufficiently simple that the integral in the last term can be done analytically, in dimensional 
regularization, so that the IR poles it generates can be canceled against those from the loop 
term. 

To build a set of universal terms that parametrize the IR singularities of any amplitude, we 
start from the observation that gauge theory amplitudes factorize in the soft limit, as follows: 




where parton j is a soft gluon, partons i, j, and k form a chain of color-space index contractions 
(we say they are color-connected), gs is the strong coupling, and the terms in parenthesis are 
called the soft eikonal factor. We here show it including mass corrections, which appear if i 
and k have non-zero rest masses, with the invariants Sab then defined as 

Sab = 2pa ■ Pb = (Pa + Pb)^ - ml - TJlf, . (44) 

The color factor, Nc, is valid for the leading-color contribution, regardless of whether the 
i and k partons are quarks or gluons. At subleading color, an additional soft-eikonal factor 
identical to the one above but with a color factor proportional to —1/Nc arises for each qq 
pair combination. This, e.g., modifies the effective color factor for qq — )- qgq from Nc to 
Nc{l — l/Nc) = 2Cf, in agreement with the color factor for quarks being Cp rather than Ca- 

Similarly, amplitudes also factorize in the collinear limit (partons i and j parallel, so 
Sij — )• 0), in which the eikonal factor above is replaced by the famous Altarelli-Parisi splitting 
kernels [34], which were already mentioned in section 2.2, in the context of PDF evolution. 
They are also the basis of conventional parton-shower models, such as those in Pythia [62]. 
We return to parton showers in section 3.2. 

Essentially, what antenna functions, CS dipoles, and the like, all do, is to combine the soft 
(eikonal) and collinear (Altarelli-Parisi) limits into one universal set of functions that achieve 
the correct limiting behavior for both soft and collinear radiation. To give an explicit example. 
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the antenna function for gluon emission from a color-connected qq pair can be derived from 
the matrix elements squared for the process ^ qq ^ qgq [63], 



\M{Z^ 



9i '^Cf 



SIK 



— + 
S 



Sjk 



collincar 



(45) 



where we have neglected mass corrections (see [64, 65] for massive expressions) and we 
recognize the universal eikonal soft factor from equation (43) in the first term. The two 
additional terms are less singular, and are required to obtain the correct collinear (Altarelli- 
Parisi) limits as Sij — ;> or Sjk — ?• 0. 

However, since the singularity structure is universal, we could equally well have used the 



process 



qq 



qgq to derive the antenna function. Our antenna function would then 



have come out as [65], 

\M{H^ qjgjqk) 
\M{H^^qiqK)\ 



Sij Sjk Sjk 
eikonal 



collincar 



2 

Sjk J SlK^ 

finite 



■^jfc _|_ ^ I _|_ 
Sii 



(46) 



where the additional term 2/sik is non-singular ("finite") over all of phase space. Thus, we 
here see an explicit example that the singularities are process-independent while the non- 
singular terms are process-dependent. Since we add and subtract the same term in equa- 
tion (42), the final answer does not depend on the choice of finite terms. We say that they 
correspond to different subtraction schemes. One standard antenna subtraction scheme, which 
uses the antenna function defined in equation (45) rather than the one in equation (46), is 
the Gehrmann-Gehrmann-de Ridder-Glover (GGG) one, given in [56]. 

If there is more than one color antenna in the Born-level process, the form of das is 
obtained as a sum over terms, each of which captures one specific soft limit and either all or 
"half" of a collinear one, depending on the specific scheme and the type of parton. 



das = ^iK-,ijk \Mf{- ■■,I,K, 



(47) 



with the sum running over all singular 3—^2 "clusterings" of the {F + l)-parton state to F 
partons. An analysis of the different ways of partitioning the collinear singularity of gluons 
among neighboring antenna is beyond the scope of this introduction, but useful discussions 
can be found in [66, 67, 68] . 



2.5 Infrared Safety 

A further requirement for being able to perform calculations within perturbative QCD is that 
the observable be infrared safe. Note: by "infrared", we here mean any limit that involves a 
low scale (i.e., any non-UV limit), without regard to whether it is collinear or soft. 

The property of infrared safety defines a special class of observables which have minimal 
sensitivity to long-distance physics, and which can be consistently computed in pQCD. An 
observable is infrared safe if: 
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1. (Safety against soft radiation): Adding any number of infinitely soft particles should not 
change the value of the observable. 

2. (Safety against collinear radiation): Splitting an existing particle up into two comoving 
particles, with arbitrary fractions z and 1 — z, respectively, of the original momentum, 
should not change the value of the observable. 

If both of these conditions are satisfied, any long-distance non-perturbative corrections will be 
suppressed by the ratio of the long-distance scale to the short-distance one to some (observable- 
dependent) power, t5^ically 



where Quv denotes a generic hard scale in the problem, and Qir ~ Aqcd ~ 0(1 GeV). 

Due to this power suppression, IR safe observables are not so sensitive to our lack of ability 
to solve the strongly coupled IR physics, unless of course we go to processes for which the 
relevant hard scale, Quy, is small (such as minimum-bias, soft jets, or small-scale jet substruc- 
ture). Even when a high scale is present, however, as in resonance decays, jet fragmentation, 
or underlying-event-type studies, infrared safety only guarantees us that infrared corrections 
are small, not that they are zero. Thus, ultimately, we run into a precision barrier even for IR 
safe observables, which only a reliable understanding of the long-distance physics itself can 
address. 

To constrain models of long-distance physics, one needs infrared sensitive observables. 
Hence it is not always the case that infrared safe observables are preferable — the purpose 
decides the tool. Instead of the suppressed corrections above, the perturbative prediction for 
such observables contains logarithms of the type already encountered in equation (36), 



which grow increasingly large as Qir/Qvy 0. As an example, consider such a fundamental 
quantity as particle multiplicities (= number of particles); in the absence of nontrivial infrared 
effects, the number of partons tends logarithmically to infinity as the IR cutoff is lowered. 
Similarly, the distinction between a charged and a neutral pion only occurs in the very last 
phase of hadronization, and hence observables that only include charged tracks, for instance, 
are always IR sensitive^"*. 

Two important categories of infrared safe observables that are widely used are event shapes 
and jet algorithms. Jet algorithms are perhaps nowhere as pedagogically described as in the 
2009 ESHEP lectures by Salam [49, Chapter 5]. Event shapes in the context of hadron colliders 
have not yet been as widely explored, but the basic phenomenology is introduced also by 
Salam and collaborators in [69], with first measurements reported by CMS and ATLAS [70, 71] 
and a proposal to use them also for the characterization of soft-QCD ("minimum-bias") events 
put forth in [72]. 

^"^This remains true in principle even if the tracks are clustered into jets, although the energy clustered in this 
way does provide a lower bound on Quv in the given event, since "charged + neutral > charged-only". 



IR Safe Observables: IR corrections cx 




(48) 



IR Sensitive Observables: IR Corrections oc 




(49) 
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Let us here merely emphasize that the real reason to prefer infrared safe jet algorithms 
over unsafe ones is not that they necessarily give very different or "better" answers in the ex- 
periment — experiments are infrared safe by definition, and the difference between infrared 
safe and unsafe algorithms may not even be visible when running the algorithm on experi- 
mental data — but that it is only possible to compute perturbative QCD predictions for the 
infrared safe ones. Any measurement performed with an infrared unsafe algorithm can only 
be compared to calculations that include a detailed hadronization model. This both limits the 
number of calculations that can be compared to and also adds an a priori unknown sensitivity 
to the details of the hadronization description, details which one would rather investigate and 
constrain separately, in the framework of more dedicated fragmentation studies. 

For LHC phenomenology, the preferred IR safe algorithm for jet reconstruction is currently 
the anti-kT one [73], with size parameter R var3dng between 0.4 and 0.7, though larger sizes 
can be motivated in certain contexts, e.g., to look for highly energetic jets and/or the boosted 
decay products of high-mass objects [74, 13]. This algorithm generates circular-looking jets, 
so subtracting off the energy believed to be associated with the underlying event (UE, see 
section 5.2) is particularly simple. 

For jet substructure, typically either the "kT" or "Cambridge/ Aachen" algorithms are used, 
see e.g. [74, 13]. The clustering measures used in these algorithms more closely mimic the 
singularity structure of QCD bremsstrahlung and they are therefore particularly well suited to 
"unravel" a tree of QCD branchings, such as a parton shower generates. The Cambridge/ Aachen 
algorithm may also be used to characterize the underlying event, see [75]. 
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3 Monte Carlo Event Generators 

In this section, we discuss the physics of Monte Carlo event generators and their mathematical 
foundations, at an introductory level. We shall attempt to convey the main ideas as clearly 
as possible without burying them in an avalanche of technical details. References to more 
detailed discussions are included where applicable. We assume the reader is already familiar 
with the contents of the preceding section on hard processes. 

The task of a Monte Carlo event generator is to calculate everything that happens in a high- 
energy collision, from the hard short-distance physics to the long wavelengths of hadronization 
and hadron decays. Obviously, this requires some compromises to be made. General-purpose 
generators like Herwig, Pythia, and Sherpa, start from low-order (LO or NLO) descriptions 
of the perturbative hard physics and then attempt to include the "most significant" corrections, 
such as higher-order matrix-element corrections and parton showers, resonance decays and 
finite-width effects, underlying event, beam remnants, hadronization, and hadron decays. 
Each of them had slightly different origins, which carries through to the emphasis placed on 
various physics aspects today: 

• Pythia. Successor to Jetset (begun in 1978). Originated in hadronization studies. 
Main feature: the Lund string fragmentation model. 

• Herwig. Successor to Earwig (begun in 1984). Originated in perturbative coherence 
studies. Main feature: angular-ordered parton showers. 

• Sherpa. Begun in 2000. Originated in studies of the matching of hard-emission matrix 
elements with parton showers. Main feature: CKKW matching. 

There is also a large number of more specialized generators, mainly for hard processes within 
and beyond the SM, including Alpgen, Calchep, Comphep, Madgraph, Whizard, and oth- 
ers, and a few that offer alternative shower models, such as Ariadne and Vincia. The most 
common interface between hard-process and parton-shower generators is the Les Houches 
Event File (LHEF) standard, defined in [76, 77] and "spoken" by most modern generator 
tools. 

Hard processes were the topic of section 2. In this section, we shall focus mainly on parton 
showers, with some brief comments on resonance decays at the end. Section 4 then concerns 
the matching of matrix elements and parton showers. Finally, models of hadronization and 
the underlying event are the topic of section 5. 

Several of the discussions below rely on material from the section on Monte Carlo Event 
Generators in the PDG Review of Particle Physics [14] and on the more comprehensive review 
by the MCnet collaboration in [78]. The latter also contains brief descriptions of the physics 
implementations of each of the main general-purpose event generators on the market, together 
with a guide on how to use (and not use) them in various connections, and a collection of 
comparisons to important experimental distributions. We highly recommend readers to obtain 
a copy of that review, as it is the most comprehensive and up-to-date review of event generators 
currently available. Another useful and pedagogical review on event generators is contained 
in the 2006 ESHEP lectures by Torbjorn Sjostrand [28], with a more recent update in [79]. 
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Relative uncertainty with n points 1-Dim d-Dim ricvai/point 



Trapezoidal Rule Xjr? Xjr?!^ 2'^ 

Simpson's Rule l/n^ 3'^ 



Monte Carlo 1/^^ 1/^ 1 



Table 2: Relative uncertainty after n evaluations, in 1 and d dimensions, for two traditional 
numerical integration methods and stochastic Monte Carlo. The last column shows the number 
of function evaluations that are required per point, in d dimensions. 

3.1 The Monte Carlo Method 

A ubiquitous problem in fundamental physics is the following: given a source located some 
distance from a detector, predict the number of counts that should be observed within the solid 
angle spanned by the detector (or within a bin of its phase-space acceptance), as a function 
of the properties of the source, the intervening medium, and the efficiency of the detector 
Essentially, the task is to compute integrals of the form 



with da a differential cross section for the process of interest. 

In particle physics, phase space has three dimensions per final-state particle (minus four 
for overall four-momentum-conservation). Thus, for problems with more than a few outgoing 
particles, the dimensionality of phase space increases rapidly. At LEP, for instance, the total 
multiplicity of neutral -I- charged hadrons (before weak decays) was typically ~ 30 particles, 
for about 86 dimensions. 

If we try to generalize the standard numerical-integration methods that are highly efficient 
in ID, we find that most of them give very slow convergence rates for higher-dimensional 
problems. For illustration, a table of convergence rates in 1 and d dimensions is given in 
table 2, comparing the Trapezoidal (2-point) rule and Simpson's (3-point) rule to random- 
number-based Monte Carlo. In ID, the 1/n^ convergence rate of the Trapezoidal rule is much 
faster than the stochastic 1 / ^/n of random-number Monte Carlo, and Simpson's rule converges 
even faster However, as we go to d dimensions, the convergence rate of the n-point rules all 
degrade with d (while the number of function evaluations required for each "point" simultane- 
ously increases). The MC convergence rate, on the other hand, remains the simple stochastic 
l/-y/n, independent of d, and each point still only requires one function evaluation. These 
are some of the main reasons that MC is the preferred numerical integration technique for 
high-dimensional problems. In addition, the random phase-space vectors it generates can be 
re-used in many ways, for instance as input to iterative solutions, to compute many different 
observables simultaneously, and/ or to hand "events" to propagation and detector-simulation 
codes. 

Therefore, virtually all numerical cross section calculations are based on Monte Carlo tech- 
niques in one form or another, the simplest being the Rambo algorithm [80] which can be 
expressed in about half a page of code and generates a flat scan over n-body phase space^^. 

^^Strictly speaking, Rambo is only truly uniform for massless particles. Its massive variant makes up for phase- 




da_ 

dn 



(50) 
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"This risk, that convergence is only given 
with a certain probability, is inherent in 
Monte Carlo calculations and is the rea- 
son why this technique was named after the 
world's most famous gambling casino. In- 
deed, the name is doubly appropriate because 
the style of gambling in the Monte Carlo 
casino, not to be confused with the noisy and 
tasteless gambling houses of Las Vegas and 
Reno, is serious and sophisticated. " 

F. James, "Monte Carlo theory and practice". 
Kept. Prog. Phys. 43 (1980) 1145 

Figure 15: Left: The casino in Monaco. Right: extract from [84] concerning the nature of 
Monte Carlo techniques. 

However, due to the infrared singularities in perturbative QCD, and due to the presence of 
short-lived resonances, the functions to be integrated, IMp+kl"^, can be highly non-uniform, 
especially for large k. This implies that we will have to be clever in the way we sample phase 
space if we want the integration to converge in any reasonable amount of time — simple 
algorithms like Rambo quickly become inefficient for k greater than a few. To address this bot- 
tleneck, the simplest step up from Rambo is to introduce generic (i.e., automated) importance- 
sampling methods, such as offered by the Vegas algorithm [81, 82]. This is still the dominant 
basic technique, although most modern codes do employ several additional refinements, such 
as several different copies of Vegas running in parallel (multi-channel integration), to further 
optimize the sampling. Alternatively, a few algorithms incorporate the singularity structure 
of QCD explicitly in their phase-space sampling, either by directly generating momenta dis- 
tributed according to the leading-order QCD singularities, in a sort of "QCD-preweighted" 
analog of Rambo, called Sarge [83], or by using all-orders Markovian parton showers to 
generate them (ViNCiA [67, 68]). 

The price of using random numbers is that we must generalize our notion of convergence. 
In calculus, we say that a sequence {A} converges to B if an n exists for which the difference 
|^i>n — ^1 < e for any e > 0. In random- number-based techniques, we cannot completely rule 
out the possibility of very pathological sequences of "dice rolls" leading to large deviations 
from the true goal, hence we are restricted to say that {A} converges to i? if an n exists for 
which the probability for |yli>n — B\ < e, for any e > 0, is greater than P, for any P G [0, 1] [84] . 
This risk, that convergence is only given with a certain probability, is the reason why Monte 
Carlo techniques were named after the famous casino in Monaco, illustrated in figure 15. 

3.2 Theoretical Basis of Parton Showers 

In section 2, we noted two conditions that had to be valid for fixed-order truncations of the 
perturbative series to be valid: firstly, the strong coupling as must be small for perturbation 
theory to be valid at all. This restricts us to the region in which all scales Qi » Aqcd- 

space biases by returning weighted momentum configurations. 
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We shall maintain this restriction in this section, i.e., we are still considering perturbative QCD. 
Secondly, however, in order to be allowed to truncate the perturbative series, we had to require 
cTfc ^ fjfc, i.e., the corrections at successive orders must become successively smaller, which 
— due to the enhancements from soft/coUinear singular (conformal) dynamics — effectively 
restricted us to consider only the phase-space region in which all jets are "hard and well- 
separated", equivalent to requiring all Qi/Qj ~ 1. In this section, we shall see how to lift this 
restriction, extending the applicability of perturbation theory into regions that include scale 
hierarchies, Qi » Qj » Aqcdj such as occur for soft jets, jet substructure, etc. 

In fact, the simultaneous restriction to all resolved scales being larger than Aqcd 
no large hierarchies is extremely severe, if taken at face value. Since we collide and observe 
hadrons (— ;> low scales) while simultaneously wishing to study short-distance physics processes 
(— )• high scales), it would appear trivial to conclude that fixed-order pQCD is not applicable to 
collider physics at all. So why do we still use it? 

The answer lies in the fact that we actually never truly perform a fixed-order calculation in 
QCD. Let us repeat the factorized formula for the cross section, equation (34), now inserting 
also a function, D, to represent the fragmentation of the final-state partons into observable 
hadrons, 

^ = E ['dx,dxjJ2 [d<^f f^/hA^^,^^hfJ/h,{xJ,^^l)^^^Df{^^o,^^j,) , (si) 

i,j f 

with O denoting the observable evaluated on the partonic final state, and O the observable 
evaluated on the hadronic final state, after fragmentation. Although the partonic cross sec- 
tion, daij^f , does represent a fixed-order calculation, the parton densities, /j/^j^ and fj/h2J 
include so-called resummations of perturbative corrections to all orders from the initial scale 
of order the mass of the proton, up to the factorization scale, fip (see section 2.2 and/or the 
TASI lectures by Sterman [30]). Note that the oft-stated mantra that the PDFs are purely 
non-perturbative functions is therefore misleading. True, they are defined as essentially non- 
perturbative functions at some very low scale, /io ~ a few GeV, but, if hf is taken large, they 
necessarily incorporate a significant amount of perturbative physics as well. On the "fixed- 
order side", all we have left to ensure in daij^f is then that there are no large hierarchies 
remaining between fxp and the QCD scales appearing in Likewise, in the final state, the 
fragmentation functions, Df, include infinite-order resummations of perturbative corrections 
all the way from fip down to some low scale, with similar caveats concerning mantras about 
their non-perturbative nature as for the PDFs. 

3.2.1 Step One: Infinite Legs 

The infinite-order resummations that are included in objects such as the PDFs and FFs in 
equation (51) (and in their parton-shower equivalents) rely on some very simple and powerful 
properties of gauge field theories that were already touched on in section 2. In particular, we 
saw in section 2.4 that we can represent all the IR limits of any NLO amplitude with a set of 
simple universal functions, based solely on knowing which partons are color-connected (i.e., 
have color-space index contractions) with one another 

The diagrams in figure 16 show the basic origin of the universal IR singularities of gauge 
theory amplitudes. On the left is shown a diagram (squared) in which an emission with 
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Figure 16: Diagrams giving rise to coUinear (left) and soft (right) singularities. 



Figure 17: Illustration of the branching phase space for qq — > qgq, with the original dipole- 
antenna oriented horizontally, the two parents sharing the transverse component of recoil, and 
the azimuthal angle (p (representing rotations of the emitted parton around the dipole axis) 
chosen such that the gluon is radiated upwards. From [67]. 

small Sij interferes with itself. In the coUinear limit, Sij — ^ 0, the propagator of the parent 
parton, /, goes on shell; the singularity of the associated propagator factor is the origin of the 
1/sij coUinear singularities. On the right is shown the interference between emission from 
parton / and emission from parton K. The resulting term has propagator singularities when 
both partons / and K go on shell, which can happen simultaneously if parton j is soft. This 
generates the 2sik/{sijSjk) soft singularity, also called the soft eikonal factor or the dipole 
factor. 

We now understand the fundamental origin of the IR singularities, why they are universal, 
and why amplitudes factorize in the soft and coUinear limits — the singularities are simply 
generated by intermediate parton propagators going on shell, which is independent of the 
nature of the hard process, and hence can be factorized from it. 

Thus, for each pair of (massless) color-connected partons / and K in F, the squared am- 
plitude for F + 1 gluon, |A^f+iPj will include a factor 



where g1 = Anas is the strong coupling, i and k represent partons / and K after the branching 
(i.e., they include possible recoil effects) and Sij is the invariant between parton i and the 
emitted parton, j. 

The branching phase space of a color dipole (i.e., a pair of partons connected by a color- 
index contraction) is illustrated in figure 17. Expressed in the branching invariants, Sij and sjk, 
the phase space has a characteristic triangular shape, imposed by the relation s = Sij + sjk + Sik 
(assuming massless partons). Sketchings of the post-branching parton momenta have been 
inserted in various places in the figure, for illustration. The soft singularity is located at the 
origin of the plot and the coUinear regions lie along the axes. 




(52) 



Antenna Function 
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The coUinear terms for aqq ^ qgq "antenna" are unambiguous and are given in section 2.4. 
Since gluons are in the adjoint representation, they carry both a color and an anticolor index 
(one corresponding to the rows and the other to the columns of the Gell-Mann matrices), 
and there is therefore some ambiguity concerning how to partition coUinear radiation among 
the two antennae they participate in. This is discussed in more detail in [67]. Differences 
are subleading, however, and for our purposes here we shall consider gluon antenna ends as 
radiating just like quark ones. The difference between quark and gluon radiation then arise 
mainly because gluons participate in two antennae, while quarks only participate in one. This 
is related to the difference between the color factors, Ca ^ 2Cf- 

The problem that plagued the fixed-order truncations in section 2 is clearly visible in equa- 
tion (52): if we integrate over the entire phase space including the region Sij — )• 0, Sjk — 0, 
we end up with a double pole. If we instead regulate the divergence by cutting off the inte- 
gration at some minimal perturbative cutoff scale fif^, we end up with a logarithm squared of 
that scale. This is a typical example of "large logarithms" being generated by the presence of 
scale hierarchies. Also note that the precise definition of ^uir is not unique. Any scale choice 
that properly isolates the singularities from the rest of phase space will do, with some typical 
choices being, for example, invariant-mass and/or transverse-momentum scales. 

Before we continue, it is worth noting that equation (52) is often rewritten in other forms 
to emphasize specific aspects of it. One such rewriting is thus to reformulate the invariants Sij 
appearing in it in terms of energies and angles, 

Sij = 2EiEj (1 - cos %) . (53) 

Rewritten in this way, the differentials can be partial-fractioned, 

^^^^^ ^ ^^^^ , ^^^^ 

This kind of rewriting enables an intuitively appealing categorization of the singularities as 
related to vanishing energies and angles, explaining why they are called soft and coUinear, 
respectively. Arguments based on this rewriting have led to important insights in QCD. For 
instance, within the framework of conventional parton showers, it was shown in a sequence 
of publications (see [85, 86] and references therein) that the destructive interference effects 
between two or more color-connected partons (coherence) can be described by using the angle 
of the emissions as the shower ordering variable. One should still keep in mind, however, 
that Lorentz non-invariant formulations come with similar caveats and warnings as do gauge 
non-invariant formulations of quantum field theory: while they can be practical to work with 
at intermediate stages of a calculation, one should be careful with any physical conclusions 
that rely explicitly on them. 

We shall therefore here restrict ourselves to a Lorentz-invariant formalism based directly 
on equation (52), pioneered by the dipole formulation of QCD cascades [63]. The coUinear 
limit is then replaced by a more general single-pole limit in which a single parton-parton in- 
variant vanishes (as, /or instance, when a pair of partons become coUinear), while the soft limit 
is replaced by one in which two (or more) invariants involving the same parton vanish simul- 
taneously (as, for instance by that parton becoming soft in a frame defined by two or more 
hard partons) . This avoids frame-dependent ambiguities from entering into the language, at 
the price of a slight reinterpretation of what is meant by coUinear and soft. 
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Figure 18: A selection of parton-shower evolution variables, represented as contours over the 
dipole phase space. Note: the right-most variable corresponds to evolution of only one of the 
parents, the one with no collinear singularity along the bottom of the plot. 



In the generator landscape, angular ordering is used by the Herwig [86] and Herwig++ [87] 
programs, and an angular veto is imposed on the virtuality-ordered evolution in Pythia 6 [88]. 
Variants of the dipole approach is used by the Ariadne [89], Sherpa [90, 91], and Vincia [92] 
programs, while the pi^-ordered showers in Pythia 6 and 8 represent a hybrid, combining 
collinear splitting kernels with dipole kinematics [93]. Phase-space contours of equal value of 
some of these choices are illustrated in figure 18. During the shower evolution, each model ef- 
fectively "sweeps" over phase space in the order implied by these contours. E.g., a p^-ordered 
dipole shower (leftmost plot in figure 18) will treat a hard-collinear branching as occurring 
"earlier" than a soft one, while a mass-ordered dipole shower (second plot) will tend to do the 
opposite. This affects the tower of virtual corrections generated by each shower model via the 
so-called Sudakov factor, discussed below. 

Independently of rewritings and philosophy, the real power of equation (52) lies in the fact 
that it is universal. Thus, for any process F, we can apply equation (52) in order to get an 
approximation for dap+i ■ We may then, for instance, take our newly obtained expression for 
F + 1 as our arbitrary process and crank equation (52) again, to obtain an approximation for 
daF+2 , and so forth. What we have here is therefore a very simple recursion relation that can 
be used to generate approximations to leading-order cross sections with arbitrary numbers of 
additional legs. The quality of this approximation is governed by how many terms besides the 
leading one shown in equation (43) are included in the game. Including all possible terms, 
the most general form for the cross section at F + n jets, restricted to the phase-space region 
above some infrared cutoff scale ^ir, has the following algebraic structure, 

4°!^ = (In^" + ln2"-i + In^"-^ + . . . + In + ^) (55) 

where we use the notation In'^ without an argument to denote generic functions of transcen- 
dentalitry A (the logarithmic function to the power A being a "typical" example of a function 
with transcendentality A appearing in cross section expressions, but also dilogarithms and 
higher logarithmic functions^^ of transcendentality > 1 should be implicitly understood to be- 
long to our notation In'^) . The last term, J", represents a rational function of transcendentality 
0. We shall also use the nomenclature singular and finite for the In'*' and T terms, respectively, 
a terminology which reflects their respective behavior in the limit —?- 0. 

^^Note: due to the theorems that allow us, for instance, to rewrite dilogarithms in different ways with logarithmic 
and lower "spillover" terms, the coefficients at each A are only well-defined up to reparametrization ambiguities 
involving the terms with transcendentality greater than A. 
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Figure 19: Coefficients of the perturbative series covered by LO + LL approximations to higher- 
multiplicity tree-level matrix elements. Green (darker) shading represents the full perturbative 
coefficient at the respective k and £. Yellow (lighter) shading represents an LL approximation 
to it. Half-shaded boxes indicate phase spaces in which we are prohibited from integrating 
over the IR singular region, as discussed in sections 2.3 and 4. 

The simplest approximation one can build on equation (55), dropping all but the leading 
In^^ term in the parenthesis, is thus the leading-transcendentality approximation. This approx- 
imation is better known as the DLA (double logarithmic approximation), since it generates the 
correct coefficient for terms which have two powers of logarithms for each power of as, while 
terms of lower transcendentalities are not guaranteed to have the correct coefficients. In so- 
called LL (leading-logarithmic) parton shower algorithms, one generally expects to reproduce 
the correct coefficients for the In^" and In^""'^ terms. In addition, several formally subleading 
improvements are normally also introduced in such algorithms (such as explicit momentum 
conservation, gluon polarization and other spin-correlation effects [94, 95, 96], higher-order 
coherence effects [85], renormalization scale choices [97], finite-width effects [98], etc), as 
a means to improve the agreement with some of the more subleading coefficients as well, if 
not in every phase-space point then at least on average. Though LL showers do not magically 
acquire NLL (next-to-leading-log) precision from such procedures, one therefore still expects 
a significantly better average performance from them than from corresponding "strict" LL ana- 
lytical resummations. A side effect of this is that it is often possible to "tune" shower algorithms 
to give better-than-nominal agreement with experimental distributions, by adjusting the pa- 
rameters controlling the treatment of subleading effects. One should remember, however, that 
there is a limit to how much can be accomplished in this way — at some point, agreement 
with one process will only come at the price of disagreement with another, and at this point 
further tuning would be meaningless. 

Applying such an iterative process on a Born-level cross section, one obtains the description 
of the full perturbative series illustrated in figure 19. The yellow (lighter) shades used here for 
k > 1 indicate that the coefficient obtained is not the exact one, but rather an approximation to 
it that only gets its leading singularities right. However, since this is still only an approximation 
to infinite-order tree-level cross sections (we have not yet included any virtual corrections), we 
cannot yet integrate this approximation over all of phase space, as illustrated by the yellow 
boxes being only half filled in figure 19; otherwise, the summed total cross section would still 
be infinite. This particular approximation would therefore still appear to be very useless indeed 
— on one hand, it is only guaranteed to get the singular terms right, but on the other, it does 
not actually allow us to integrate over the singular region. In order to obtain a truly all-orders 
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Figure 20: Coefficients of the perturbative series covered by LO + LL calculations, impos- 
ing unitarity order by order for each n = k + i. Green (darker) shading represents the full 
perturbative coefficient at the respective k and £. Yellow (lighter) shading represents an LL 
approximation to it. 



calculation, the constraint of unitarity must also be explicitly imposed, which furnishes an 
approximation to all-orders loop corrections as well. Let us therefore emphasize that figure 19 
is included for pedagogical purposes only; all resummation calculations, whether analytical 
or parton-shower based, include virtual corrections as well and consequently yield finite total 
cross sections, as will now be described. 

3.2.2 Step Two: Infinite Loops 

Order-by-order unitarity, such as used in the KLN theorem, implies that the singularities caused 
by integration over unresolved radiation in the tree-level matrix elements must be canceled, 
order by order, by equal but opposite-sign singularities in the virtual corrections at the same 
order. That is, from equation (52), we immediately know that the 1-loop correction to dap 
must contain a term, 

2Re[MPM?*] D -glNc 

that cancels the divergence coming from equation (52) itself. Further, since this is universally 
true, we may apply equation (56) again to get an approximation to the corrections generated 
by equation (52) at the next order and so on. By adding such terms explicitly, order by order, 
we may now bootstrap our way around the entire perturbative series, using equation (52) to 
move horizontally and equation (56) to move along diagonals of constant n = k + I. Since 
real-virtual cancellations are now explicitly restored, we may finally extend the integrations 
over all of phase space, resulting in the picture shown in figure 20. 

The picture shown in figure 20, not the one in figure 19, corresponds to what is actually 
done in resummation calculations, both of the analytic and parton-shower types^^. Physically, 
there is a significant and intuitive meaning to the imposition of unitarity, as follows. 

Take a jet algorithm, with some measure of jet resolution, Q, and apply it to an arbitrary 
sample of events, say dijets. At a very crude resolution scale, corresponding to a high value 

^^In the way these calculations are formulated in practice, they in fact rely on one additional property, called 
exponentiation, that allows us to move along straight vertical lines in the loops-and-legs diagrams. However, since 
the two different directions furnished by equations (52) and (56) are already sufficient to move freely in the full 
2D coefficient space, we shall use exponentiation without extensively justifying it here. 



M^p' / — ^i_^]L I — _|_ less singular terms ) , (56) 
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for Q, you find that everything is clustered back to a dijet configuration, and the 2-jet cross 
section is equal to the total inclusive cross section, 

O"tot = 0"F;incl ■ (57) 

At finer resolutions, decreasing Q, you see that some events that w^ere previously classified as 
2-jet events contain additional, lower-scale jets, that you can now resolve, and hence those 
events now migrate to the 3-jet bin, while the total inclusive cross section of course remains 
unchanged, 

0"tot = CrF;mcI = <7i?;excl(Q) + 0'F+l;incl (Q) , (58) 

where "incl" and "excl" stands for inclusive and exclusive cross sections^^, respectively, and the 
Q-dependence in the two terms on the right-hand side must cancel so that the total inclusive 
cross section is independent of Q. Later, some 3-jet events now migrate further, to 4 and higher 
jets, while still more 2-jet events migrate into the 3-jet bin, etc. For arbitrary n and Q, we have 

n-l 

^F+n;mc\{Q) = Cr_F;incl " ^ 0-F+m;cxcliQ) ■ (59) 

m=0 

This equation expresses the trivial fact that the cross section for n or more jets can be computed 
as the total inclusive cross section for F minus a sum over the cross sections for F + exactly 
m jets including all m < n. On the theoretical side, it is these negative terms which must 
be included in the calculation, for each order n = A; + ^, to restore unitarity. Physically, they 
express that, at a given scale Q, each event will be classified as having either 0, 1, 2, or 
whatever jets. Or, equivalently for each event we gain in the 3-jet bin as Q is lowered, we 
must loose one event in the 2-jet one; the negative contribution to the 2-jet bin is exactly 
minus the integral of the positive contribution to the 3-jet one, and so on. We may perceive 
of this detailed balance as an evolution of the event structure with Q, for each event, which is 
effectively what is done in parton-shower algorithms, to which we shall return in Section 3.3. 



3.3 Perturbation Theory with Markov Chains 



Consider again the Born-level cross section for an arbitrary hard process, F, differentially in 
an arbitrary infrared-safe observable O, as obtained from equation (35): 



dof 
dO 



Born 



(60) 



where the integration runs over the full final-state on-shell phase space of F (this expression 
and those below would also apply to hadron collisions were we to include integrations over 
the parton distribution functions in the initial state), and the 6 function projects out a 1- 
dimensional slice defined by O evaluated on the set of final-state momenta which we denote 

To make the connection to parton showers, we insert an operator, S, that acts on the 
Born-level final state before the observable is evaluated, i.e., 

dap 



dO 



S 



d<^F \Mf\^S{<!>F,0) . 



(61) 



inclusive = F plus anything. F exclusive = F and only F. Thus, aF;inci = X^k^o ^F+fc; 



+ fc;cxcl 
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Formally, this operator — the evolution operator — will be responsible for generating all 
(real and virtual) higher-order corrections to the Born-level expression. The measurement 5 
function appearing explicitly in equation (60) is now implicit in S. 

Algorithmically parton showers cast S as an iterative Markov (i.e., history-independent) 
chain, with an evolution parameter, Qe, that formally represents the factorization scale of 
the event, below which all structure is summed over inclusively. Depending on the partic- 
ular choice of shower algorithm, Q e may be defined as a parton virtuality (virtuality-order 
showers), as a transverse-momentum scale (p^-ordered showers), or as a combination of en- 
ergies times angles (angular ordering). Regardless of the specific form of Qe, the evolution 
parameter will go towards zero as the Markov chain develops, and the event structure will 
become more and more exclusively resolved. A transition from a perturbative evolution to a 
non-perturbative one can also be inserted, when the evolution reaches an appropriate scale, 
typically around 1 GeV. This scale, called the hadronization scale, thus represents the lowest 
perturbative scale that can appear in the calculations, with all perturbative corrections below 
it summed over inclusively. 

Working out the precise form that S must have in order to give the correct expansions 
discussed in section 3.2 takes a bit of algebra, and is beyond the scope we aim to cover in these 
lectures. Heuristically the procedure is as follows. We noted that the singularity structure of 
QCD is universal and that at least its first few terms are known to us. We also saw that we 
could iterate that singularity structure, using universality and unitarity, thereby bootstrapping 
our way around the entire perturbative series. This was illustrated by figure 20 in section 3.2. 

Skipping intermediate steps, the form of the all-orders pure-shower Markov chain, for the 
evolution of an event between two scales Qei > Qe2> is, 

S{<^p, Qei,Qe2, O) = A(«>^, Qei.Qe2) 6 {O - 0{^e)) 



F + exclusive above Q e2 

Sr{^F+i) Qei, Qf+i) S{^f+i,Qf+i, Qe2, O) 



E2 '^^E 



F + \ inclusive above Q e2 

(62) 



with the so-called Sudakov factor, 

A{^e,Qei,Qe2) = exp 



r JQe2 ^^E 



(63) 



defining the probability that there is no evolution (i.e., no emissions) between the scales Qei 
and Qe2, according to the radiation functions Sr to which we shall return below. The term on 
the first line of equation (62) thus represents all events that did not evolve as the resolution 
scale was lowered from Qei to Qe2, while the second line contains a sum and phase-space in- 
tegral over those events that did evolve — including the insertion of S{^p^^) representing the 
possible further evolution of the event and completing the iterative definition of the Markov 
chain. 

The factor d<I>^^^ /d*!*^ defines the chosen phase space factorization. Our favorite is the 
so-called dipole-antenna factorization, whose principal virtue is that it is the simplest Lorentz 
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invariant factorization which is simultaneously exact over all of phase space while only involv- 
ing on-shell momenta. For completeness, its form is 

d^F+i m . . d<A 1 ,,,, 
^ = d^ = "^'^^ ^^^^ 2^ 16^ ' ^^4^ 

which involves just one color-anticolor pair for each r, with invariant mass squared Sr = 
{Pa + Pi + Pb)"^- Other choices, such as purely coUinear ones (only exact in the coUinear 
limit or involving explicitly off-shell momenta), more global ones involving all partons in the 
event (more complicated, in our opinion), or less global ones with a single parton playing 
the dominant role as emitter, are also possible, again depending on the specific algorithm 
considered. 

The radiation functions Sr obviously play a crucial role in these equations, driving the 
emission probabilities. For example, if Sr — )• 0, then A — )• exp(O) = 1 and all events stay in the 
top line of equation (62). Thus, in regions of phase space where Sr is small, there is little or 
no evolution. Conversely, for Sr — ^ oo, we have A — ^ 0, implying that all events evolve. One 
possible choice for the radiation functions Sr was implicit in equation (52), in which we took 
them to include only the leading (double) singularities, with r representing color-anticolor 
pairs. In general, the shower may exponentiate the entire set of universal singular terms, or 
only a subset of them (for example, the terms leading in the number of colors Nc), which is 
why we here let the explicit form of Sr be unspecified. Suffice it to say that in traditional parton 
showers, Sr would simply be the DGLAP splitting kernels (see, e.g., [4]), while they would be 
so-called dipole or antenna radiation functions in the various dipole-based approaches to QCD 
(see, e.g., [63, 50, 56, 91, 67, 68]). 

The procedure for how to technically "construct" a shower algorithm of this kind, using 
random numbers to generate scales distributed according to equation (62), is described more 
fully in [67], using a notation that closely parallels the one used here. The procedure is 
also described at a more technical level in the review [78], though using a slightly different 
notation. Finally, a pedagogical introduction to Monte Carlo methods in general can be found 
in [84]. 



3.4 Decays of Unstable Particles 

In most BSM processes and some SM ones, an important aspect of the event simulation is how 
decays of short-lived particles, such as top quarks, EW and Higgs bosons, and new BSM reso- 
nances, are handled. We here briefly summarize the spectrum of possibilities, but emphasize 
that there is no universal standard. Users are advised to check whether the treatment of a 
given code is adequate for the physics study at hand. 

The appearance of an unstable resonance as a physical particle at some intermediate stage 
of the event generation implies that its production and decay processes are treated as being 
factorized. This is called the narrow width approximation and is valid up to corrections of 
order F/mo, with T the width and mo the pole mass of the particle. States whose widths are a 
substantial fraction of their mass should not be treated in this way, but rather as intrinsically 
virtual internal propagator lines. 

For states treated as physical particles, two aspects are relevant: the mass distribution of 
the decaying particle itself and the distributions of its decay products. For the mass distribu- 
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tion, the simplest is to use a 5 function at mo. The next level up, typically used in general- 
purpose Monte Carlos, is to use a Breit-Wigner distribution (relativistic or non-relativistic), 
which formally resums higher-order virtual corrections to the mass distribution. Note, how- 
ever, that this still only generates an improved picture for moderate fluctuations away from mo. 
Similarly to above, particles that are significantly off-shell (in units of F) should not be treated 
as resonant, but rather as internal off-shell propagator lines. In most Monte Carlo codes, some 
further refinements are also included, for instance by letting F be a function of m ("running 
widths") and by limiting the magnitude of the allowed fluctuations away from mo. See also 
[99] for an elaborate discussion of the Higgs boson lineshape. 

For the distributions of the decay products, the simplest treatment is again to assign them 
their respective mo values, with a uniform (i.e., isotropic, or "flat") phase-space distribution. 
A more sophisticated treatment distributes the decay products according to the differential 
decay matrix elements, capturing at least the internal d3mamics and helicity structure of the 
decay process, including EPR-like correlations. Further refinements include polarizations of 
the external states [94, 95, 96] (see also [100, 101, 102] for phenomenological studies) and 
assigning the decay products their own Breit-Wigner distributions, the latter of which opens 
the possibility to include also intrinsically off-shell decay channels, like H — ^ WW*. Please 
refer to the physics manual of the code you are using and/ or make simple cross checks like 
plotting the distribution of phase-space invariants it produces. 

During subsequent showering of the decay products, most parton-shower models will pre- 
serve the total invariant mass of each resonance-decay system separately, so as not to skew the 
original resonance shape. 
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Figure 21: The double-counting problem caused by naively adding cross sections involving 
matrix elements with different numbers of legs. 



4 Matching at LO and NLO 

The essential problem that leads to matrix-element/parton-shower matching can be illustrated 
in a very simple w^ay. Assume we have computed the LO cross section for some process, F, 
and that we have added an LL shower to it, as in the left-hand pane of figure 21. We know 
that this only gives us an LL description of F + 1. We now wish to improve this from LL to LO 
by adding the actual LO matrix element for F + 1. Since we also want to be able to hadronize 
these events, etc, we again add an LL shower off them. However, since the matrix element for 
F + 1 is divergent, we must restrict it to cover only the phase-space region with at least one 
hard resolved jet, illustrated by the half-shaded boxes in the middle pane of figure 21. 

Adding these two samples, however, we end up counting the LL terms of the inclusive cross 
section for F + 1 twice, since we are now getting them once from the shower off F and once 
from the matrix element for F + 1, illustrated by the dark shaded (red) areas of the right- 
hand pane of figure 21. This double-counting problem would grow worse if we attempted to 
add more matrix elements, with more legs. The cause is very simple. Each such calculation 
corresponds to an inclusive cross section, and hence naive addition would give 

0"tot = CO;incl + <7l;mcl = O'Ojexcl + 2cri;incl • (65) 

Recall the definition of inclusive and exclusive cross sections, equation (58): F inclusive = F 
plus anything. F exclusive = F and only F. Thus, aF;mci = J2kLo(^F+k;exc\- 

Instead, we must match the coefficients calculated by the two parts of the full calculation 
— showers and matrix elements — more systematically, for each order in perturbation theory, 
so that the nesting of inclusive and exclusive cross sections is respected without overcounting. 

Given a parton shower and a matrix-element generator, there are fundamentally three 
different ways in which we can consider matching the two [67]: slicing, subtraction, and 
unitarity. The following subsections will briefly introduce each of these. 

4.1 Slicing 

The most commonly encountered matching type is currently based on separating (slicing) 
phase space into two regions, one of which is supposed to be mainly described by hard matrix 
elements and the other of which is supposed to be described by the shower. This type of ap- 
proach was first used in Herwig [103], to include matrix- element corrections for one emission 
beyond the basic hard process [104, 105]. This is illustrated in figure 22. The method has 
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Figure 22: Herwig's original matching scheme [104, 105], in which the dead zone of the 
Herwig shower was used as an effective "matching scale" for one emission beyond a basic 
hard process. 



since been generalized by several independent groups to include arbitrary numbers of addi- 
tional legs, the most well-known of these being the CKKW [106], CKKW-L [107, 108], and 
MLM [109, 110] approaches. 

Effectively, the shower approximation is set to zero above some scale, either due to the 
presence of explicit dead zones in the shower, as in Herwig, or by vetoing any emissions 
above a certain matching scale, as in the (L)-CKKW and MLM approaches. The empty part of 
phase space can then be filled by separate events generated according to higher-multiplicity 
tree-level matrix elements (MEs). In the (L)-CKKW and MLM schemes, this process can be 
iterated to include arbitrary numbers of additional hard legs (the practical limit being around 
3 or 4, due to computational complexity) . 

In order to match smoothly with the shower calculation, the higher-multiplicity matrix 
elements must be associated with Sudakov form factors (representing the virtual corrections 
that would have been generated if a shower had produced the same phase-space configura- 
tion), and their as factors must be chosen so that, at least at the matching scale, they become 
identical to the choices made on the shower side. The CKKW and MLM approaches do this by 
constructing "fake parton-shower histories" for the higher-multiplicity matrix elements. By ap- 
plying a sequential jet clustering algorithm, a tree-like branching structure can be created that 
at least has the same dominant structure as that of a parton shower Given the fake shower 
tree, as factors can be chosen for each vertex with argument as{p±) and Sudakov factors can 
be computed for each internal line in the tree. In the CKKW method, these Sudakov factors 
are estimated analytically, while the MLM and CKKW-L methods compute them numerically, 
from the actual shower evolution. 

Thus, the matched result is identical to the matrix element (ME) in the region above the 
matching scale, modulo higher-order (Sudakov and as) corrections. We may sketch this as 

ME corrections 
Matched (above matching scale) = Exact x {l + 0{as)) , (66) 

where the "shower-corrections" include the approximate Sudakov factors and reweighting 
factors applied to the matrix elements in order to obtain a smooth transition to the shower- 
dominated region. 

Below the matching scale, the small difference between the matrix elements and the 
shower approximation can be dropped (since their leading singularities are identical and this 
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Figure 23: Slicing, with up to two additional emissions beyond the basic process. The showers 
off F and F + 1 are set to zero above a specific "matching scale". (The number of coefficients 
shown was reduced a bit in these plots to make them fit in one row.) 



region by definition includes no hard jets), yielding the pure shower answer in that region, 

shower correction 

X N ^ N 

Matched (below matching scale) = Approximate + (Exact — Approximate) 

= Approximate + non-singular 
— ;> Approximate . (67) 

This type of strategy is illustrated in figure 23. 

As emphasized above, since this strategy is discontinuous across phase space, a main point 
here is to ensure that the behavior across the matching scale be as smooth as possible. CKKW 
showed [106] that it is possible to remove any dependence on the matching scale through 
NLL precision by careful choices of all ingredients in the matching; technical details of the 
implementation (affecting the 0{as) terms in eq. (66)) are important, and the dependence 
on the unphysical matching scale may be larger than NLL unless the implementation matches 
the theoretical algorithm precisely [107, 108, 111]. Furthermore, since the Sudakov factors 
are generally computed using showers (MLM, L-CKKW) or a shower-like formalism (CKKW), 
while the real corrections are computed using matrix elements, care must be taken not to (re- 
)introduce differences that could break the detailed real-virtual balance that ensures unitarity 
among the singular parts, see e.g., [112]. 

It is advisable not to choose the matching scale too low. This is again essentially due 
to the approximate scale invariance of QCD imploring us to write the matching scale as a 
ratio, rather than as an absolute number. If one uses a very low matching scale, the higher- 
multiplicity matrix elements will already be quite singular, leading to very large LO cross 
sections before matching. After matching, these large cross sections are tamed by the Sudakov 
factors produced by the matching scheme, and hence the final cross sections may still look 
reasonable. But the higher-multiplicity matrix elements in general contain subleading singu- 
larity structures, beyond those accounted for by the shower, and hence the delicate line of 
detailed balance that ensures unitarity has most assuredly been overstepped. We therefore 
recommend not to take the matching scale lower than about an order of magnitude below the 
characteristic scale of the hard process. 

One should also be aware that all strategies of this type are quite computing intensive. 
This is basically due to the fact that a separate phase-space generator is required for each of 
the n-parton correction terms, with each such sample a priori consisting of weighted events 
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Figure 24: Mc@nlo. In the middle pane, cyan boxes denote non-singular correction terms, 
while the egg-colored ones denote showers off such corrections, which cannot lead to double- 
counting at the LL level. 



such that a separate unweighting step (often with quite low efficiency) is needed before an 
unweighted sample can be produced. 

4.2 Subtraction 

Another way of matching two calculations is by subtracting one from the other and correcting 
by the difference, schematically 

shower correction 
Matched = Approximate + (Exact — Approximate) . (68) 

This looks very much like the structure of a subtraction-based NLO fixed-order calculation, 
section 2.4, in which the shower approximation here plays the role of subtraction terms, and 
indeed this is what is used in strategies like Mc@nlo [113, 114, 115], illustrated in figure 24. 
In this type of approach, however, negative-weight events will generally occur, for instance in 
phase-space points where the approximation is larger than the exact answer. 

Negative weights are not in principle an insurmountable problem. Histograms can be filled 
with each event counted according to its weight, as usual. However, negative weights do affect 
efficiency. Imagine a worst-case scenario in which 1000 positive-weight events have been 
generated, along with 999 negative-weight ones (assuming each event weight has the same 
absolute value, the closest one can get to an unweighted sample in the presence of negative 
weights) . The statistical precision of the MC answer would be equivalent to one event, for 
2000 generated, i.e., a big loss in convergence rate. 

In practice, generators like MC@NLO "only" produce around 10% or less events with nega- 
tive weights, so the convergence rate should not be severely affected for ordinary applications. 
Nevertheless, the problem of negative weights motivated the development of the so-called 
POWHEG approach [116], illustrated in figure 25, which is constructed specifically to prevent 
negative-weight events from occurring and simultaneously to be more independent of which 
parton-shower algorithm it is used with. In the Powheg method, one effectively modifies the 
real-emission probability for the first emission to agree with the F + 1 matrix element (this is 
covered under unitarity below). One is then left with a purely virtual correction, which will 
typically be positive, at least for processes for which the NLO cross section is larger than the 
LO one. 
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Figure 25: Powheg. In the middle pane, cyan boxes denote non-singular correction terms, 
while the egg-colored ones denote showers off such corrections, which cannot lead to double- 
counting at the LL level. 
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Figure 26: Menlops. Note that each of the Powheg and CKKW samples are composed of 
separate sub-samples, as illustrated in figures 23 and 25. 



The advantage of these methods is obviously that NLO corrections to the Born level can 
be systematically incorporated. However, a systematic way of extending this strategy beyond 
the first additional emission is not available, save for combining them with a slicing-based 
strategy for the additional legs, as in Menlops [117], illustrated in figure 26. These issues 
are, however, no more severe than in ordinary fixed-order NLO approaches, and hence they 
are not viewed as disadvantages if the point of reference is an NLO computation. 



4.3 Unitarity 

The oldest, and in my view most attractive, approach [88, 118] consists of working out the 
shower approximation to a given fixed order, and correcting the shower splitting functions at 
that order by a multiplicative factor given by the ratio of the matrix element to the shower 
approximation, phase-space point by phase-space point. We may sketch this as 

, correction 
shower , ^ , 



Exact 



Matched = Approximate x -, . (69) 

Approximate 

When these correction factors are inserted back into the shower evolution, they guarantee 
that the shower evolution off n — 1 partons correctly reproduces the n-parton matrix elements, 
without the need to generate a separate n-parton sample. That is, the shower approximation 
is essentially used as a pre-weighted (stratified) all-orders phase-space generator, on which a 
more exact answer can subsequently be imprinted order by order in perturbation theory. Since 
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Figure 27: Pythia (left) and Vincia (right). Unitarity-based. Only one event sample is pro- 
duced by each of these methods, and hence no sub-components are shown. 



the shower is already optimized for exactly the kind of singular structures that occur in QCD, 
very fast computational speeds can therefore be obtained with this method [68]. 

In the original approach [88, 118], used by Pythia [119, 62], this was only worked out 
for one additional emission beyond the basic hard process. In Powheg [116, 120], it was 
extended to include also virtual corrections to the Born-level matrbc element. Finally, in Vin- 
cia [67, 68], it has been extended to include arbitrary numbers of emissions at tree level, 
though that method has so far only been applied to final-state showers. 

An illustration of the perturbative coefficients that can be included in each of these ap- 
proaches is given in figure 27, as usual with green (darker shaded) boxes representing exact 
coefficients and yellow (light shaded) boxes representing logarithmic approximations. 

Finally, two more properties unique to this method deserve mention. Firstly, since the 
corrections modify the actual shower evolution kernels, the corrections are automatically re- 
summed in the Sudakov exponential, which should improve the logarithmic precision once 
> 2 is included, and secondly, since the shower is unitary, an initially unweighted sam- 
ple of (n — l)-parton configurations remains unweighted, with no need for a separate event- 
unweighting or event-rejection step. 
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5 Hadronization and Soft Hadron-Hadron Physics 

We here give a very brief overview of the main aspects of soft QCD that are relevant for 
hadron-hadron collisions, such as hadronization, minimum-bias and soft-inclusive physics, 
and the so-called underlying event. This will be kept at a pedestrian level and is largely based 
on the reviews in [14, 78, 121]. 

In the context of event generators, hadronization denotes the process by which a set of 
colored partons (after showering) is transformed into a set of color-singlet primary hadrons, 
which may then subsequently decay further. This non-perturbative transition takes place at 
the hadronization scale Qhadj which by construction is identical to the infrared cutoff of the 
parton shower. In the absence of a first-principles solution to the relevant dynamics, event 
generators use QCD-inspired phenomenological models to describe this transition. 

Essentially, the problem that a hadronization model should address can be stated as fol- 
lows: given a set of partons resolved at a scale of Qhad ~ 1 GeV, we need a "mapping" from 
this set onto a set of on-shell color-singlet (i.e., confined) hadronic states. MC models do this 
in three steps: 

1 . Map the partonic system onto a continuum of high-mass hadronic states (called "strings" 
or "clusters"). 

2. Iteratively map strings/ clusters onto discrete set of primary hadrons (string breaks / 
cluster splittings / cluster decays) . 

3. Sequential decays into secondary hadrons (e.g., p — )• vrvr, A — )• rnr, vr*^ — > 77, ...). 

The physics governing this mapping is non-perturbative and hence cannot (yet?) be solved 
from first principles. However, we do have some knowledge of the properties that such a 
solution must have. For instance, Poincare invariance, unitarity and causality are all concepts 
that apply beyond perturbation theory. In addition, lattice QCD provides us a means of making 
explicit quantitative studies in a genuinely non-perturbative setting (albeit only of certain 
questions, see section 1). 

An important result in "quenched" lattice QCD^^ is that the potential of the color-dipole 
field between a charge and an anticharge appears to grow linearly with the separation of the 
charges, at distances greater than about a femtometer This is known as "linear confinement", 
and it forms the starting point for the string model of hadronization, discussed below in sec- 
tion 5.1. Alternatively, a property of perturbative QCD called "preconfinement" [122] is the 
basis of the cluster model of hadronization, described in [78, 14]. 

In the generator landscape, Pythia uses string fragmentation, while Herwig and Sherpa 
use cluster fragmentation. It should be emphasized that the so-called parton level that can 
be obtained by switching off hadronization in an MC generator, is not a universal concept, 
since each model defines the hadronization scale differently. E.g., the hadronization scale can 
be defined by a cutoff in invariant mass, transverse momentum, or some other quantity, with 
different tunes using different values for the cutoff. Comparisons to distributions at this level 
(i.e., with hadronization switched off) may therefore be used to provide an idea of the overall 
impact of hadronization corrections within a given model, but should be avoided in the context 
of physical observables. 

^'Quenched QCD implies no "dynamical" quarks, i.e., no g ~^ qq splittings allowed. 
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Figure 28: Illustration of the transition between a Coulomb potential at short distances to the 
string-like one of equation (70) at large qq separations. 



We use the term "soft hadron-hadron physics" to comprise all scattering processes for which 
a hard, perturbative scale is not required to be present^°. This includes elastic, diffractive, 
minimum-bias, and pile-up processes, as well as the physics contributing to the so-called un- 
derlying event. We give a brief introduction to such processes in section 5.2. 

We round off with a discussion of the data constraints that enter in the tuning of Monte 
Carlo models in section 5.3, and give an outline of a procedure that could be followed in a 
realistic set-up. 

5.1 The String Model of Hadronization 

Starting from early concepts developed by Artru and Mennessier [124], several hadronization 
models based on strings were proposed in the late 1970'ies and early 80'ies. Of these, the 
most widely used today is the so-called Lund model, implemented in the Pythia code. We 
shall therefore concentrate on that particular model here, though many of the overall concepts 
would be shared by any string-inspired method. (A more extended discussion can be found in 
the very complete and pedagogical review of the Lund model by Andersson [125].) 

Consider the production of a qq pair from vacuum, for instance in the process e+e" — ^ 
'y*/Z —?- qq ^ hadrons. As the quarks move apart, linear confinement implies that a potential 

V{r) = Kr (70) 

is asymptotically reached for large distances, r. (At short distances, there is a Coulomb term 
proportional to 1/r as well, but this is neglected in the Lund model.) This potential describes 
a string with tension (energy per unit length) k. The physical picture is that of a color flux 
tube being stretched between the q and the q, figure 28. From hadron mass spectroscopy the 
string tension k, is known to be 

K ~ IGeV/fm ~ 0.2 GeV^. (71) 

A straightforward Lorentz-invariant description of this object is provided by the massless rela- 
tivistic string in 1-1-1 dimensions, with no transverse degrees of freedom. The mathematical, 
one-dimensional string can be thought of as parameterizing the position of the axis of a cylin- 
drically symmetric flux tube. (Note that the expression "massless" is somewhat of a misnomer, 
since k effectively corresponds to a "mass density" along the string.) 

^"Note, however, that while a hard scale is not required to be present, it is not explicitly required to be absent 
either. Thus, both diffractive, minimum-bias, pile-up and underlying-event processes will have tails towards high- 
Pi physics as well. For example, even tt pair production can be viewed as a tail of minimum-bias interactions, and 
there is a tail of diffractive processes in which hard dijets can be produced diffractively (see, e.g., [123]). 
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leftover string, 
further breaks 




Figure 29: a) Illustration of string breaking by quark pair creation in the string field, b) 
Illustration of the algorithmic choice to process the fragmentation from the outside-in, splitting 
off a single on-shell hadron in each step. 



As the q and q move apart, their kinetic energy is gradually converted to potential en- 
ergy, stored in the growing string spanned between them. In the "quenched" approximation, 
in which g ^ qq splittings are not allowed, this process would continue until the endpoint 
quarks have lost all their momentum, at which point they would reverse direction and be ac- 
celerated by the now shrinking string. In the real world, quark-antiquark fluctuations inside 
the string field can make the transition to become real particles by absorbing energy from 
the string, thereby screening the original endpoint charges from each other and breaking the 
string into two separate color-singlet pieces, {qq) — > {qq') + {q'q), illustrated in figure 29a. This 
process then continues until only ordinary hadrons remain. (We will give more details on the 
individual string breaks below.) 

More complicated multi-parton topologies including gluons are treated by representing 
gluons as transverse "kinks". Thus soft gluons effectively "build up" a transverse structure in 
the originally one-dimensional object, with infinitely soft ones absorbed into the string without 
leading to modifications. For strings with finite-energy kinks, the space-time evolution is then 
slightly more involved [125], and modifications to the fragmentation model to handle stepping 
across gluon corners have to be included, but the main point is that there are no separate 
free parameters for gluon jets. Differences with respect to quark fragmentation arise simply 
because quarks are only connected to a single string piece, while gluons have one on either 
side, increasing the energy loss per unit (invariant) time from a gluon to the string by a factor 
of 2 relative to quarks, which can be compared to the ratio of color Casimirs Ca/Cf = 2.25. 

Since the string breaks are causally disconnected (as can easily be realized from space-time 
diagrams [125]), they do not have to be considered in any specific time-ordered sequence. In 
the Lund model, the string breaks are instead generated starting with the leading ("outer- 
most") hadrons, containing the endpoint quarks, and iterating inwards towards the center of 
the string, alternating randomly between fragmentation off the left- and right-hand sides, re- 
spectively, figure 29 b. This has the advantage that a single on-shell hadron can be split off 
in each step, making it straightforward to ensure that only states consistent with the known 
spectrum of hadron resonances are produced, as will be discussed below. 

The details of the individual string breaks are not known from first principles. The Lund 
model invokes the idea of quantum mechanical tunneling, which leads to a Gaussian suppres- 
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sion of the energies and masses imparted to the produced quarks, 

PToh{m^,p]_ ) cx exp I 1 exp I — 1 , (72) 




where ruq is the mass of the produced quark and p± is the transverse momentum imparted to 
it by the breakup process (the antiquark obviously has the same mass and opposite p^). 

Due to the factorization of the p±and m dependence implied by equation (72), the spectrum 
of produced quarks in this model is independent of the quark flavor, with a universal average 
value of 

(pig) = = k/tt {250MeVf . (73) 

Bear in mind that "transverse" is here defined with respect to the string axis. Thus, the p±m a 
frame where the string is moving is modified by a Lorentz boost factor Also bear in mind that 
(T^ is here a purely non-perturbative parameter In a Monte Carlo model with a fixed shower 
cutoff Qhadj the effective amount of "non-perturbative" p^may be larger than this, due to 
effects of additional unresolved soft-gluon radiation below Qhad- In principle, the magnitude 
of this additional component should scale with the cutoff, but in practice it is up to the user 
to enforce this by retuning (see section 5.3) the effective a parameter when changing the 
hadronization scale. Since hadrons receive pj_ contributions from two breakups, one on either 
side, their average transverse momentum squared will be twice as large, 

(Plh) = 2^' . (74) 

The mass suppression implied by equation (72) is less straightforward to interpret. Since 
quark masses are notoriously difficult to define for light quarks, the value of the strangeness 
suppression must effectively be extracted from experimental measurements, e.g., of the K/tt 
ratio, with a resulting suppression of roughly s/u ~ s/d ~ 0.2 - 0.3. Inserting even compara- 
tively low values for the charm quark mass in equation (72), however, one obtains a relative 
suppression of charm of the order of 10^^^. Heavy quarks can therefore safely be considered 
to be produced only in the perturbative stages and not by the soft fragmentation. 

Baryon production can be incorporated in the same basic picture [126], by allowing string 
breaks to occur also by the production of pairs of so-called diquarks, loosely bound states of 
two quarks in an overall 3 representation (e.g., red -I- blue = antigreen). Again, the rela- 
tive rate of diquark-to-quark production is not known a priori and must be extracted from 
experimental measurements, e.g., of the p/ir ratio. More advanced scenarios for baryon pro- 
duction have also been proposed, in particular the so-called popcorn model [127, 128], which 
is normally used in addition to the diquark picture and then acts to decrease the correla- 
tions among neighboring baryon-antibaryon pairs by allowing mesons to be formed inbetween 
them. Within the Pythia framework, a fragmentation model including explicit string junc- 
tions [129] has so far only been applied to baryon-number-violating new-physics processes 
and to the description of beam remnants (and then acts to increase baryon stopping [130]). 

This brings us to the next step of the algorithm: assignment of the produced quarks 
within hadron multiplets. The fragmenting quark (antiquark) may combine with the anti- 
quark (quark) from a newly created breakup to produce either a vector or a pseudoscalar 
meson, or, if diquarks are involved, either a spin- 1/2 or spin-3/2 baryon. Unfortunately, the 
string model is entirely unpredictive in this respect, and this is therefore the sector that con- 
tains the largest amount of free parameters. From spin counting alone, one would expect the 
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Figure 30: Illustration of the Lund symmetric fragmentation function, correctly normalized to 
unity. Left: variation of the a parameter, from 0.1 (blue) to 0.9 (red), with fixed 6 = 1 GeV~^ 
and m_L = 1 GeV. Right: variation of the b parameter, from 0.5 (red) to 2 (blue) GeV'^, with 
fixed a = 0.5 and m± = 1 GeV. 



ratio V/ S of vectors to pseudoscalars to be 3, but in practice this is only approximately true 
for B*/B. For lighter flavors, the difference in phase space caused by the V-S mass splittings 
implies a suppression of vector production. Thus, for D* /D, the effective ratio is already re- 
duced to about ~ 1.0 - 2.0, while for K* /K and p/n, extracted values range from 0.3 - 1.0. 
Recall, as always, that these are production ratios of primary hadrons, hence feed-down (from 
secondary decays of heavier hadrons) complicates the extraction of these parameters from ex- 
perimental data, in particular for the lighter hadron species. The production of higher meson 
resonances is assumed to be low in a string framework^^. For diquarks, separate parameters 
control the relative rates of spin-1 diquarks vs. spin-0 ones and, likewise, have to extracted 
from data, with resulting values of order {qq)i/{qq)o ~ 0.075 - 0.15. 

With and now fixed, the final step is to select the fraction, z, of the fragmenting end- 
point quark's longitudinal momentum that is carried by the created hadron. In this respect, the 
string picture is substantially more predictive than for the flavor selection. Firstly, the require- 
ment that the fragmentation be independent of the sequence in which breakups are considered 
(causality) imposes a "left- right symmetry" on the possible form of the fragmentation function, 
f{z), with the solution 

/(.)ocl(l-.)'^exp(-^^M±i?L)^ ^ (75) 

which is known as the Lund symmetric fragmentation function (normalized to unit integral). 
The a and b parameters, illustrated in figure 30, are the only free parameters of the fragmen- 
tation function, though a may in principle be flavor-dependent. Note that the explicit mass 
dependence in f{z) implies a harder fragmentation function for heavier hadrons (in the rest 
frame of the string) . 

As a by-product, the probability distribution in invariant time r of q'q' breakup vertices, 
or equivalently F = (kt)^, is also obtained, with dP/dV oc r"exp(— ftF) implying an area law 

^^The four L — 1 multiplets are implemented in Pythia, but are disabled by default, largely because several 
states are poorly known and thus may result in a worse overall description when included. 
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Figure 31: Illustration of the iterative selection of flavors and momenta in the Lund string 
fragmentation model. 



for the color flux [131], and the average breakup time lying along a hyperbola of constant 
invariant time tq ~ 10~^^s [125]. 

The iterative selection of flavors, p±, and z values is illustrated in figure 31. A parton 
produced in a hard process at some high scale Quv emerges from the parton shower, at the 
hadronization scale Qir, with 3-momentum p = {p±o,P+), where the "+" on the third com- 
ponent denotes "light-cone" momentum, p± = E ± pz- Next, an adjacent pair from the 
vacuum is created, with relative transverse momenta ±p±i. The fragmenting quark combines 
with the (i from the breakup to form a 7r+, which carries off a fraction zi of the total lightcone 
momentum The next hadron carries off a fraction Z2 of the remaining momentum, etc. 

For massive endpoints (e.g., c and b quarks, or hypothetical hadronizing new-physics par- 
ticles, generally called i?-hadrons), which do not move along straight lightcone sections, the 
exponential suppression with string area leads to modifications of the form [132], f{z) — > 
with ruQ the mass of the heavy quark. Strictly speaking, this is the only fragmen- 
tation function that is consistent with causality in the string model, though a few alternative 
forms are typically provided as well. 

Note, however, that the term fragmentation function in the context of non-perturbative 
hadronization models is used to denote only the corrections originating from scales below the 
infrared cutoff scale of the parton shower. That is, the fragmentation functions introduced here 
are defined at an intrinsically low scale of order Q ~ 1 GeV. It would therefore be inconsis- 
tent to compare them directly to those that are used in fixed-order / analytical-resummation 
contexts, which are typically defined at a factorization scale of order the scale of the hard 
process. 

5.2 Soft Hadron-Hadron Processes 
5.2.1 Elastic and Inelastic Scattering 

Elastic hadron-hadron scattering consists of all reactions of the t3^e 



AipA)BipB) ^ Aip'A)B{p's) 



(76) 
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where A and B are hadrons with momenta pA and ps, respectively. Specifically, the only 
exchanged quantity is momentum; all quantum numbers and masses remain unaltered, and 
no new particles are produced. Inelastic scattering covers everything else, i.e., 

AB^X ^AB , i77) 

where X ^ AB signifies that one or more quantum numbers are changed, and/or more 
particles are produced. The total hadron-hadron cross section can thus be written as a sum of 
these two physically distinguishable components, 

O-tot(s) = CTelis) + O-inel(s) , (78) 

where s = {pa +Pb)'^ is the beam-beam center-of-mass energy squared. 
5.2.2 Diffractive Scattering 

If A and/or B are not elementary, the inelastic final states may be further divided into "diffrac- 
tive" and "non-diffractive" topologies. This is a qualitative classification, usually based on 
whether the final state looks like the decay of an excitation of the beam particles (diffrac- 
tive^^), or not (non-diffractive), or upon the presence of a large rapidity gap somewhere in 
the final state which would separate such excitations. 

Given that an event has been labeled as diffractive, either within the context of a theoreti- 
cal model, or by a final-state observable, we may distinguish between three different classes of 
diffractive topologies, which it is possible to distinguish between physically, at least in princi- 
ple. In double-diffractive (DD) events, both of the beam particles are diffractively excited and 
hence none of them survive the collision intact. In single-diffractive (SD) events, only one of 
the beam particles gets excited and the other survives intact. The last diffractive topology is 
central diffraction (CD), in which both of the beam particles survive intact, leaving an excited 
system in the central region between them. (This latter topology includes "central exclusive 
production" where a single particle is produced in the central region.) That is, 

O-inol(s) = CrSD(s) + O-Dd(s) + crCD(s) + CJnd(s) , (79) 

where "ND" (non-diffractive, here understood not to include elastic scattering) contains no 
gaps in the event consistent with the chosen definition of diffraction. Further, each of the 
diffractively excited systems in the events labeled SD, DD, and CD, respectively, may in princi- 
ple consist of several subsystems with gaps between them. Eq. (79) may thus be defined to be 
exact, within a specific definition of diffraction, even in the presence of multi-gap events. Note, 
however, that different theoretical models almost always use different (model-dependent) def- 
initions of diffraction, and therefore the individual components in one model are in general 
not directly comparable to those of another. It is therefore important that data be presented 
at the level of physical observables if unambiguous conclusions are to be drawn from them. 

^^An example of a process that would be labeled as diffractive would be if one the protons is excited to a A"*" 
which then decays back to + vr", without anything else happening in the event. In general, a whole tower of 
possible diffractive excitations are available, which in the continuum limit can be described by a mass spectrum 
falling roughly as AKP jhP. 
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5.2.3 Minimum Bias 

The term minimum-bias (MB) is an experimental term, used to define a certain class of events 
that are selected with the minimum possible trigger bias, to ensure they are as inclusive as pos- 
sible^-^. In theoretical contexts, the term "minimum-bias" is often used with a slightly different 
meaning; to denote specific (classes of) inclusive soft-QCD subprocesses in a given model. 
Since these two usages are not exactly identical, in these lectures we have chosen to reserve 
the term "minimum bias" to pertain strictly to definitions of experimental measurements, and 
instead use the term "soft inclusive" physics as a generic descriptor for the class of processes 
which generally dominate the various experimental "minimum-bias" measurements in theo- 
retical models. This parallels the terminology used in the review [78], from which most of 
the discussion here has been adapted. See equation (79) above for a compact overview of the 
t5^es of physical processes that contribute to minimum-bias data samples. For a more detailed 
description of Monte Carlo models of this physics, see [78]. 

5.2.4 Underljdng Event and Multiple Parton Interactions 

In events containing a hard parton-parton interaction, the underlying event (UE) can be roughly 
conceived of as the difference between QCD with and without including the remnants of the 
original beam hadrons. Without such "beam remnants", only the hard interaction itself, and 
its associated parton showers and hadronization, would contribute to the observed particle 
production. In reality, after the partons that participate in the hard interaction have been taken 
out, the remnants still contain whatever is left of the incoming beam hadrons, including also 
a partonic substructure, which leads to the possibility of multiple parton interactions (MPI). 
Useful reviews of MPI-based MC models can be found in [78, 121]. Analytical models are 
mostly formulated only for double parton scattering, see e.g., [133, 134, 135, 136]. 

Due to the simple fact that the remnants are not empty, an underlying event will always 
be there — but how much additional energy does it deposit in a given measurement region? 
A quantification of this can be obtained, for instance, by comparing measurements of the UE 
to the average activity in minimum-bias events at the same ^/s. Interestingly, it turns out 
that the underlying event is much more active, with larger fluctuations, than the average 
MB event. This is called the jet pedestal effect (hard jets sit on top of a higher-than-average 
"pedestal" of underlying activity), and is interpreted as follows: when two hadrons collide at 
non-zero impact parameter, high-p^ interactions can only take place inside the overlapping 
region. Imposing a hard trigger therefore statistically biases the event sample toward more 
central collisions, which will also have more underlying activity. 

For hard processes at the LHC at 7 TeV, the transverse energy, Et, in the UE is about 1.5 
GeV per unit AR area, though with large event-to-event fluctuations of order ±1 GeV [137]. 
Thus, using the typical sizes for LHC jets quoted in section 2.5, and multiplying by ttR"^, the Et 

^^A typical min-bias trigger would thus be the requirement of at least one measured particle in a given rapidity 
region, so that all events which produce at least one observable particle would be included, which must, indeed, 
be considered the minimal possible bias. In principle, everything is a subset of minimum-bias, including both hard 
and soft processes. However, compared to the total minimum-bias cross section, the fraction that is made up of 
hard processes is only a very small tail. Since only a tiny fraction of the total minimum-bias rate can normally 
be stored, the minimum-bias sample would give quite poor statistics if used for hard physics studies. Instead, 
separate dedicated hard-process triggers are typically included in addition to the minimum-bias one, in order to 
ensure maximal statistics also for hard physics processes. 
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originating from the UE, in a cone with radius 0.4, can be estimated to be Et\je{R = 0.4) ~ 
0.75±0.5 GeV, while the Et in cones with radii 0.7 and 1.0 would be Etue{R = 0.7) ~ 2.3±1.5 
GeV and Etve{R = 1-0) ~ 4.7 it 3 GeV, respectively. See [75] for a discussion on the use of jet 
algorithms in the characterization of the UE (including illustrations with simple analytical toy 
models) and these lecture notes by R. Field [138] for more discussion on the level of the UE 
at the Tevatron and at LHC. 

5.3 Tuning 

The main virtue of general-purpose Monte Carlo event generators is their ability to provide a 
complete and fully differential picture of collider final states, down to the level of individual 
particles. This allows them to be used as detailed — albeit approximate — theoretical ref- 
erences for measurements performed at accelerators like the LHC, against which models of 
both known and 'new' physics can be tested. As has been emphasized in these lectures, the 
achievable accuracy depends both on the inclusiveness of the chosen observable and on the 
sophistication of the simulation itself. An important driver for the latter is obviously the de- 
velopment of improved theoretical models, e.g., by including matching to higher-order matrix 
elements, more accurate resummations, or better non-perturbative models, as discussed in the 
previous sections; but it also depends crucially on the available constraints on the remaining 
free parameters of the model. Using existing data to constrain these is referred to as generator 
tuning. 

Although Monte Carlo models may appear to have a bewildering array of independently 
adjustable parameters, it is worth keeping at the front of one's mind that most of these pa- 
rameters only control relatively small (exclusive) details of the event generation. The majority 
of the (inclusive) physics is determined by only a few, very important ones, such as, e.g., the 
value of the strong coupling, in the perturbative domain, and the form of the fragmentation 
function for massless partons, in the non-perturbative one. 

Armed with a good understanding of the underlying model, an expert would therefore 
normally take a highly factorized approach to constraining the parameters, first constraining 
the perturbative ones and thereafter the non-perturbative ones, each ordered in a measure 
of their relative significance to the overall modeling. This factorization, and carefully chosen 
experimental distributions corresponding to each step, allows one to concentrate on just a 
few parameters and distributions at a time, reducing the full parameter space to manageable- 
sized chunks. Still, each step will often involve more than one single parameter, and non- 
factorizable corrections still imply that changes made in subsequent steps can change the 
agreement obtained in previous ones by a non-negligible amount, requiring additional itera- 
tions from the beginning to properly tune the entire generator framework. 

Recent years have seen the emergence of automated tools that attempt to reduce the 
amount of both computer and manpower required for this task, for instance by making full 
generator runs only for a limited set of parameter points, and then interpolating between these 
to obtain approximations to what the true generator result would have been for any interme- 
diate parameter point, as, e.g., in the Professor tool [139, 140]. Automating the human expert 
input is of course more difficult. In the tools currently on the market, this question is addressed 
by a combination of input solicited from the generator authors (e.g., which parameters and 
ranges to consider, which observables constitute a complete set, etc) and the elaborate con- 
struction of non-trivial weighting functions that determine how much weight is assigned to 
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each individual bin and to each distribution. The field is still burgeoning, however, and future 
sophistications are to be expected. Nevertheless, at this point the overall quality of the tunes 
obtained with automated methods appear to at least be competitive with the manual ones. 
Recent examples of tunes including uncertainty variations can be found in [141, 112, 142]. 

A sketch of a reasonably complete tuning procedure, without going into details about the 
parameters that control each of these sectors in individual Monte Carlo models, would be the 
following: 

1) Keep in mind that inabilities of models to describe data is a vital part of the feedback 
cycle between theory and experiment. Also keep in mind that perturbation theory at LOxLL 
is doing very well if it gets within 10% of a given IR safe measurement. An agreement of 5% 
should be considered the absolute sanity limit, beyond which it does not make any sense to 
tune further The advent of NLO Monte Carlos may reduce these numbers slightly, but only 
for quantities for which one expects NLO precision to hold, see section 4. However, the sanity 
limit should be taken to be at least twice as large for quantities governed by non-perturbative 
physics. For some quantities, e.g., ones for which the underlying modeling is known to be poor, 
an order-of-magnitude agreement or worse may have to be accepted. Attempting to force 
Monte Carlo models to describe data far outside their domains of validity must be expected to 
produce similar side effects as attempting to turn a Fiat into a Ferrari merely by cranking up 
the engine revolutions. 

2) Final-State Radiation and Hadronization: mainly using LEP and other e+e" collider 
data. On the IR safe side, there are event shapes and jet observables, the latter including rates, 
resolutions, masses, shapes, and jet-jet correlations. On the IR sensitive side, special attention 
should be paid to the high-z tail of the fragmentation spectra, where a single hadron carries a 
large fraction of an entire jet's momentum, since this is the tail that is most likely to give "fake 
jets". Depending on the focus of the tuning, attention should also be paid to identified-particle 
rates and ratios, and to fragmentation in events containing heavy quarks and/or gluon jets. 
Usually, more weight is given to those particles that are most copiously produced, though this 
again depends on the focus. Finally, particle-particle correlations and baryon production are 
typically some of the least well constrained components of the overall modeling. The scaling 
properties of IR safe vs. IR sensitive contributions can be tested by comparing data at several 
different e+e^ collider energies. 

3) Initial-State Radiation, and so-called "Primordial^^ kx": here, one would in prin- 
ciple like to use data from DIS reactions, which are less complicated to interpret than full 
hadron-hadron collisions. However, due to difficulties in translating between the ep and pp 
environments, this is normally not what is done in practice. Instead, the main constrain- 
ing distribution is the dilepton distribution in Drell-Yan events in hadron-hadron collisions. 
For any observables containing explicit jets, be aware that the underlying event can produce 
small horizontal shifts in jet p± distributions, which may in turn result in seemingly larger- 
than-expected vertical changes if the distributions are falling sharply. Also note that the ISR 
evolution is sensitive to the choice of PDFs, with caveats as discussed in section 2.2. 

4) Initial-Final Connections: e.g., radiation from color lines connected to the initial state 
and jet broadening in hadron collider environments. This is one of the most poorly controlled 
parts of most MC models. Keep in mind that it is not directly constrained by pure final-state 

^"^Primordial kr'- an additional soft p± component that is injected on top of the p± generated by the initial-state 
shower itself, see [78, Section 7.1]. 
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observables, such as LEP fragmentation, nor by pure initial-state observables, such as the Drell- 
Yan p_i_ spectrum, which is why we list it as a separate item here. In principle, DIS would again 
be a prime territory for placing constraints on this aspect at least for quark jets, but in practice 
more often inclusive-jet and other multi-jet processes (such as W/Z+ jets) in hadron colliders 
are used. The modeling of this aspect can have important effects on specific observables, a 
recent example being the tt forward-backward asymmetry at the Tevatron [143]. 

5) Underlying Event: Good constraints on the overall level of the underlying event can be 
obtained by counting the summed transverse energy (more IR safe) and/ or particle multiplici- 
ties and average transverse momenta (more IR sensitive) in regions transverse to a hard trigger 
jet (more IR safe) or particle (more IR sensitive), see e.g. [138]. Constraints on the fluctua- 
tions of the underlying event are also important, and can be obtained, e.g., by comparing to 
measurements of the RMS of such distributions [137]. Again, note that the UE is sensitive to 
the choice of PDFs [144]. 

6) Color (Re-) Connections and other Final-State Interactions: By Final-State Interac- 
tions, we intend a broad spectrum of possible collective effects that may be included to a 
greater or lesser extent in various models. These effects include: Bose-Einstein correlations 
(see, e.g., [145]), rescattering (see, e.g., [146]), color reconnections / string interactions (see, 
e.g., [147, 148, 149]), hydrodynamics (see, e.g., [150]), etc. As a rule, these effects are soft 
and/ or non-perturbative and hence should not modify hard IR safe observables appreciably. 
They can, however, have drastic effects on IR sensitive ones, such as particle multiplicities, mo- 
mentum distributions, and correlations, wherefore useful constraints are typically furnished by 
measurements of spectra and correlations as functions of quantities believed to serve as indi- 
cators of the strength of these phenomena (such as event multiplicity), and/ or by collective- 
flow-type measurements. Finally, if the model includes a universal description of underlying 
event and soft-inclusive QCD, as many MPI-based models do, then minimum-bias data can 
also be used as a control sample, though one must then be careful either to address diffractive 
contributions properly or to include only data samples that minimize their impact. 

7) Beam Remnants: Constraints on beam remnant fragmentation (see, e.g., [130]) are 
most easily obtained in the forward region, but, e.g., the amount of baryon transport from the 
remnant to a given rapidity region can also be used to probe how much the color structure of 
the remnant was effectively disturbed, with more baryon transport indicating a larger amount 
of "beam baryon blowup". 

We round off by emphasizing that comparisons of specific models and tunes to data can be 
useful both as immediate tests of commonly used models, and to illustrate the current amount 
of theoretical uncertainty surrounding a particular distribution. Independently of how well 
the models fit the data, such comparisons also provide a set of well-defined theoretical refer- 
ence curves that serve as useful guidelines for future studies. However, the conclusions that 
can be drawn from comparisons of individual tunes of specific models on single distributions 
are necessarily limited. In order to obtain more general conclusions, a strategy for a more 
coherent and over-arching look at both the data and the models was recently proposed in 
[144]. Specifically, rather than performing one global tune to all the data, as is usually done, 
a more systematic check on the validity of the underl)dng physics model can be obtained by 
instead performing several independent optimizations of the model parameters for a range 
of different phase-space windows and/or collider environments. In regions in which consis- 
tent parameter sets are obtained, with predictions that are acceptably close to the data, the 
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underlying model can then be considered as interpolating well, i.e., it is universal. If not, 
a breakdown in the ability of the model ability to span different physical regimes has been 
identified, and can be addressed, with the nature of the deviations giving clues as to the na- 
ture of the breakdown. With the advent of automated tools making it easier to run several 
optimizations without much additional computing overhead, such systematic studies are now 
becoming feasible, with a first example given in [144]. 
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